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f- « ) Abstract. Motivated by numerical simulations showing the emergence of either pe- 

£SJ \ riodic or irregular patterns, wc explore a mechanism of pattern formation arising in 

; i ■ the processes described by a system of a single reaction-diffusion equation coupled with 

Oh| ordinary differential equations. We focus on a basic model of early canccrogenesis pro- 

<T^ ■ posed by Marciniak-Czochra and Kimmcl [Comput. Math. Methods Med. 7 (2006), 189- 

213], [Math. Models Methods Appl. Sci. 17 (2007), suppl., 1693-1719], but the theory 
we develop applies to a wider class of pattern formation models with an autocatalytic 
non-diffusing component. The model exhibits diffusion-driven instability (Turing- type 
instability). However, we prove that all Turing-type patterns, i.e., regular stationary 
solutions, are unstable in the Lyapunov sense. Furthermore, we show existence of dis- 
continuous stationary solutions, which are also unstable. 
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1. Introduction 

> '. 

. In this paper, we explore a mechanism of pattern formation arising in processes de- 

scribed by a system of a single react ion- diffusion equation coupled with ordinary differ- 
ential equations. Such systems of equations arise from modeling of interactions between 
cellular processes and diffusing growth factors. 

A rigorous derivation, using methods of asymptotic analysis (homogenization), of the 
macroscopic reaction-diffusion models describing the interplay between the nonhomoge- 
neous cellular dynamics and the signaling molecules diffusing in the intercellular space has 
been recently published in [8]. It was shown that receptor-ligand binding processes can 
^ ! be modeled by reaction-diffusion equations coupled with ordinary differential equations 
in the case when all membrane processes are homogeneous within the membrane, which 
seems to be the case in most of processes. More precisely, such receptor- based models can 
be represented by the following initial-boundary value problem 

u t = f(u,v), 

v t = DAv + g(u, v) in Q, 
d n v = on dn, 

u(x,0) = u Q (x), 
v(x,0) = v (x), 

where u and v are vectors of variables, D is a diagonal matrix with positive coefficients 
on the diagonal, the symbol d n denotes the normal derivative (no- flux condition), and Q 
is a bounded region. It was shown that, if homogeneity of the processes on the membrane 
does not hold, equations with additional integral terms are obtained, see [8]. 
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One of possible mechanisms of pattern formation in such models is the diffusion-driven 
instability (DDI), called also the Turing-type instability. Let us recall that the diffusion- 
driven instability arises in a reaction-diffusion system, when there exists a spatially ho- 
mogeneous solution, which is asymptotically stable with respect to spatially homogeneous 
perturbations, but unstable to spatially inhomogeneous perturbation. Models with DDI 
describe the process of destabilization of a stationary spatially homogeneous state and 
evolution of spatially heterogeneous structures, which converge to a spatially heteroge- 
neous steady state. The systems with DDI were prevalent in the modeling literature since 
the seminal paper of Turing [15] and have provided explanations of pattern formation in 
a variety of biological systems, see e.g. [10] and references therein. 

The existing qualitative theory of such systems is mainly focused on stability conditions 
for homogeneous steady states. Moreover, the majority of theoretical studies focus on the 
analysis of the non-degenerated reaction-diffusion systems, i. e. with a strictly positive dif- 
fusion coefficient in each equation. However, in many biological applications, it is relevant 
to consider receptor-based systems involving ordinary differential equations. An interest- 
ing class of such systems consists of only a single reaction-diffusion equation coupled with 
a system of ordinary differential equations. The existence and stability of spatially hetero- 
geneous patterns arising in models exhibiting diffusion-driven instability, but consisting of 
only one reaction-diffusion equation, is a mathematically interesting issue. Such models 
are very different from classical Turing-type models and the spatial structure of the pat- 
tern emerging from the destabilization of the spatially homogeneous steady state cannot 
be concluded based on linear stability analysis [10]. Asymptotic analysis of systems with 
such a dispersion relation seems to be an open problem. 

In particular, such systems arise in the modeling of the growth of a spatially-distributed 
cell population, which proliferation is controlled by endogenous or exogenous growth fac- 
tors diffusing in the extracellular medium and binding to cell surface as proposed by 
the first author and Kimmel in the series of recent papers [5, 6, 7]. The generic model 
has the following mathematical form of mixed react ion- diffusion and ordinary differential 
equations 

(1.1) u t = (a p (u,v) - d c ) u, 

(1.2) v t = a(u)w — dv — d b v, 

(1.3) w t = DAw — d g w — a(u)w + dv + k(u). 

Here, unknown functions u(x,t), v(x,t), w(x,t) describe the densities of cells, free and 
bound growth factor molecules, respectively, distributed over a certain bounded domain. 
In system (1.1)-(1.3), the proliferation rate a p (u, v), the growth factor binding rate a(u), 
and the production of the free growth factor by cells n{u) are given functions. We use 
the subscripts in the parameters d c , <4, d g in order to emphasise that they denote the 
degradation factors of cells, bound molecules and growth molecules, respectively. 

System (1.1)- (1.3) is usually supplemented with the homogeneous Neumann (zero flux) 
boundary conditions for the function w = w(x, t) and with nonnegative initial conditions. 

Based on preliminary mathematical analysis and using numerical simulations, the au- 
thors of [5, 6, 7] found conditions for the existence of a positive spatially homogeneous 
steady state exhibiting diffusion-driven instability. In numerical simulations instability of 
the constant steady state leads to the emergence of growth patterns concentrated around 
discrete points along the spatial coordinate, which take the mathematical form of spike- 
type spatially inhomogeneous solutions. This multifocality is as expected from the field 
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theory of carcinogenesis. Numerical simulations showed qualitatively new patterns of be- 
havior of solutions, including, in some cases, a strong dependence of the emerging pattern 
on initial conditions and quasi-stability followed by rapid growth of solutions. 

The main goal of this work is to develop mathematical theory that allows us to study 
properties of solutions to this kind of initial-boundary value problems. In this work, for 
simplicity of the exposition, we assume that cells occupy the interval x G [0, L] and we 
consider a particular form of functions the next section. Our results, however, 

hold true in a much more general case what we have systematically emphasized in our 
construction of patterns (see Section 5) and in the proof of their instability (Section 6). 

The paper is organized as follows. In Section 2, main results of this paper are for- 
mulated. Section 3 is devoted to the existence of solutions. Section 4 provides prelimi- 
nary results on the large time behavior of solutions and, in particular, a criterion on the 
pointwise extinction of solutions. In Section 5, nonhomogeneous stationary solutions are 
characterized and, in Section 6, instability of all stationary solutions is shown using linear 
stability analysis. The paper is supplemented by Appendix devoted to the analysis of the 
corresponding kinetics system (system of ordinary differential equations). 

Notation. The usual norm of the Lebesgue space L p (0, 1) is denoted by || • || p for any 
p G [l,oo] and W k,p (0,l) is the corresponding Sobolev space. The constants (always 
independent of x and t) will be denoted by the same letter C, even if they may vary from 
line to line. Sometimes, we write, e.g., C = C(p,q,r, ...) when we want to emphasize the 
dependence of C on parameters p,q,r, .... 



2. Results and coments 

2.1. Statement of the problem. In this work, we consider system (1.1)- (1.3) on the 
bounded interval x G [0,L]. Moreover, in our mathematical analysis, we assume that the 
proliferation rate has the Hill function form 

, . v/u av 

a p (u, v) = a- 



1 + v/u u + v 

with a given constant a > 0. We consider the quadratic growth factor binding rate 
a(u) = m 2 , which follows from conditions for the diffusion driven instability, see [5] and 
the next subsection for more details. Moreover, we assume that the production of the free 
growth factor is constant, k(u) = k,q > 0. This assumption is imposed for simplicity of 
the exposition; note that a more general case is also interesting from the modeling point 
of view, see [6]. Finally, we introduce the diffusion coefficient I/7, which is a composite 
parameter including the diffusion constant D > and the scaling parameter L > 0, 
namely, we set 7 = L 2 /D. 

To summarize, we study the following system of three ordinary/partial differential 
equations 

(2.1) Ut= {£~^~ d <) u for x G [0,1], t > 0, 

(2.2) v t = -d b v + u 2 w - dv for x G [0,1], t > 0, 

(2.3) w t = —w xx — dgiu — u 2 w + dv + for x G (0, 1), t > 
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supplemented with the homogeneous Neumann (zero flux) boundary conditions for the 
function w = w(x,t) 

(2.4) w x (0, t) = w x (l, t) = for all t > 
and with nonnegative initial conditions 

(2.5) u(x, 0) = Uq(x), v(x,0) — v (x), w(x, 0) = w (x). 

Here, the letters a, d c , d&, d g , d, 7, kq denote positive constants. 

In Section 3, we show that the initial-boundary value problem (2.1)-(2.5) has a unique 
and global-in-time solution for a large class of nonnegative initial conditions. Such results 
on the global-in-time existence, the regularity of solutions, and their positivity for all 
t > are rather standard for reaction-diffusion equations with non-zero diffusion in each 
equation, see e.g. [11] and the references therein. A more careful analysis is required in 
the case of the ODE-PDE system, hence, we state these results for the completeness of 
the exposition. 

2.2. Diffusion-driven instability. We are interested in systems with the diffusion- 
driven instability. As it was stated in [7, Prop. 3.1], a generic system of two ordinary 
differential equations coupled with a reaction-diffusion equation, and such that an < for 
i = 1, 2, 3 and 012021 > 0, exhibits the diffusion-driven instability if there exists a positive, 
spatially constant steady state, for which the following conditions are satisfied 



(2.6) -tr (A) > 0, 

(2.7) -tr (A) det(Ai) + det(A) > 0, 

i<j 

(2.8) -det(A)>0, 

(2.9) -det(A 2 ) > 0, 

where A = (an) . . , „ „ is the Jacobian matrix of the system without diffusion linearized 



around this constant positive equilibrium and Aij is a submatrix of A consisting of the 
i-th and j-th column and z-th and j-th row. 

Conditions (2.6)-(2.8) are necessary for the stability of the steady state in the absence of 
diffusion. Inequality (2.9) is a sufficient and necessary condition for destabilization of this 
steady state in every system with all an < for i — 1, 2, 3 and 012021 > (for the proof 
see [5]). These conditions guarantee that the model (1.1)-(1.3) exhibits diffusion- driven 
instability if the function a(u) evaluated at the steady state u satisfies a(u) — ua'(u) < 0, 
what, in particularly, always holds for a(u) = u 2 in model (2.1)-(2.3). 

In the particular case of system (2.1)-(2.3), conditions (2.6)-(2.9) hold true for the 
constant steady state u_, w_) defined in (2.12) (see the next subsection) and for the 
matrix A stated in (2.31) with W = Note that condition (2.9) leads to 2 eigenvalues 
of the matrix A12 of an opposite sign. The positive eigenvalue of ^4i 2 plays the crucial 
role in the proof of instability of all steady states of system (2.1)-(2.3), see Subsection 2.5. 

2.3. Preliminary properties of solutions. We begin our study of qualitative prop- 
erties of solutions to problem (2.1)-(2.5) by considering x-independent solutions which 
satisfy the corresponding kinetic system of the three ordinary differential equations (see 
(A.1)-(A.3), below). Let us briefly summarize our results on the kinetic system that we 
prove in Appendix A. 
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For every constant initial condition w(0) > 0, v(0) > 0, and w(0) > 0, the problem 
(2.1)-(2.5) has a unique x-independent global-in-time solution (u(t),v(t),w(t)). 
This solution stays, for t > 0, in a bounded set, see Proposition A.l. Here, if 
u(to) = v(to) = for certain to > 0, the right-hand side of equation (2.1) is 
satisfied in the limit sense, namely, when v \ and u \ 0. 
The constant vector 



(2.10) (t2o,Vo,w ) = (°' '^) 



is a trivial steady state of system (2.1)-(2.3). This is an asymptotically stable solu- 
tion not only of the kinetic system (A.1)-(A.3), but also of the react ion- diffusion 
equations (2.1)-(2.4), see Corollary 2.4 and Proposition 4.1. Moreover, in The- 
orems A.2-A.5 of Appendix, we describe convergence rates of solutions to the 
kinetic system towards the trivial steady states. 
Assume that a > d c and 

d 2 c (d b + d) 



(2.11) «5 > 6, where 6 = Ad g d 



b" 



d, 



\2 



Then, system (2.1)-(2.3) has constant positive stationary solutions (u±,v±,w±), 
where 

_ k ± \JkI - 6 _ d 2 c (d b + d) 1 _ a-d c _ 

2.12) w± = J- , v± = - — -j- — , u± = — - — v±. 

2a g [a — d c y w± a c 

• The stationary solution (u_,v_,w_) of the kinetic system (A.1)-(A.3) is asymp- 
totically stable, namely, the linearization matrix of this system at the steady state 
(m_, u>_) has all eigenvalues with negative real parts. On the other hand, 
the steady state (u + ,v + ,w + ) is an unstable solution of (A.1)-(A.3), because its 
linearization matrix has one positive eigenvalue, see Theorem A. 9 and Corollary 
A. 10, for more details. 
From now on, we consider solutions of problem (2.1)-(2.5) that are not necessarily 
space-inhomogeneous. Unlike in the case of the kinetic system (see Proposition A.l in 
Appendix), in view of numerical simulations, we do not expect that such a solution stays 
in a certain invariant region as it is discussed in the monograph by Smoller [14, Ch. 14, 
§B] . However, in the following theorem, we prove that integrals of all nonnegative solutions 
of (2.1)-(2.5) enter, as t — > oo, into a certain invariant set, which is independent of initial 
conditions. Notice that the mass of a solution is controlled for alii > as in the survey 
article [11], were systems with non-zero diffusion in each equation were discussed. 

Theorem 2.1. Assume that (u,v,w) is a nonnegative global-in-time solution of problem 
(2.1)-(2.5) corresponding to a bounded initial condition (u q ,v ,Wq) G -^ 1 (0, 1) x L 1 (0, 1) x 
L°°(0, 1). Denote fi = min{<i s , d b } > 0. Then, the following estimates hold 



K aKo 
H ' d c \i 



(2.13) limsup / u(x, t) dx < min 

t— >OD JO 

(2.14) limsup / v(x, t) dx < — , 

t-too Jo A 4 

(2.15) limsup ||iu(t)||oo < K o I 

f->°° \ fid 

with a numeric constant C > independent of problem (2.1) -(2.5). 




1/2 

'9 
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This theorem is proved at the end of Section 3. Notice that it holds also true for k = 0, 
hence, in this particular case, each nonnegative solution of problem (2.1)-(2.5) satisfies 
(IK*)||i, |K*)||i, |Kt)||oo) -»• (0,0,0) as t -)• oo. 

Next, we discuss the stability of the trivial steady state (2.10) as a solution of the 
reaction-diffusion equations (2.1)-(2.4). First of all, we have to emphasize that, under 
the assumption a < d c , each nonnegative solution of problem (2.1)-(2.5) converges expo- 
nentially towards (0,0, Ko/d g ), see Proposition 4.1 below. Hence, in this work, we have 
always to assume that a > d c to observe a nontrivial large time behavior of solutions. 

It is clear from equation (2.1) that if Uq(x) = for some x G [0, 1] than u(x,t) = for 
all t > and, by equation (2.2), we have v(x,t) — > as t — > 00. Furthermore, below in 
Proposition 4.2, we show that, for each x G [0, 1], the following two conditions 

(2.16) \im u(x, t) = and limv(x, t) = 

t— >oo t—too 

are equivalent. The following result plays a fundamental role in understanding of the 
pattern formation described by problem (2.1)-(2.5) (see Remark 2.14, below) and says 
that relations (2.16) hold true if Uq(x) and v$(x) are not too large. 

Theorem 2.2. Assume that a > d c and other coefficients in system (2.1)-(2.5) are positive 
and arbitrary. Fix x G [0, 1] and assume that there exist constants K w > and M > 
satisfying 

d c\ 2 „ , , , A ( d^ 2 



(2.17) MK W I 1 + — <(d b + d) - 

\ a ) \ a 

and such that for allt > 

(2.18) < w(x,t) < K w and < u (x) < M, < v (x) < ( ^ J M. 

Then, (u(x,t),v(x,t)) — > (0,0) as t — > 00. This is the uniform convergence for all x G 
[0, 1], for which inequalities (2.18) are satisfied. 

Remark 2.3. It follows from Theorem 2.1 that ^(t)^ is a bounded function of t > 0, 
hence, one can always find a constant K w > with the property stated in Theorem 2.2. 
In fact, this constant can be chosen explicitly, see Remark 3.5 at the end of Section 3. 

We conclude this subsection by a corollary of Theorem 2.2 on the asymptotic stability 
of the trivial steady state. 

Corollary 2.4. Let a > d c and other coefficients in system (2.1)-(2.5) be positive. Assume 
that there exist constants K w > and M > satisfying (2.17) and (2.18) for all x G [0, 1]. 
Then (u(x,t),v(x,t),w(x,t)) — > (0, 0, n /d g ) as t — > 00 uniformly inx. 

Here, it has to be emphasized that system (2.1)-(2.5) cannot be linearized around the 
solution (0, 0, K /d c ), because of the singularity at u = v = in equation (2.1). Hence, 
to prove the convergence of a solution towards the trivial steady state a direct method 
(based on the ordinary differential equations (2.1)-(2.2)) was invented, see the proof of 
Theorem 2.2 in Section 4. 
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2.4. Existence and nonexistence of patterns. Here, we describe stationary solutions 
of (2.1)-(2.4), namely, functions (U (x) , V (x) , W (x)) that satisfy the system 

(2.19) (j*L- d \u = ^ 



u + v 

(2.20) -d b V + U 2 W - dV = 0, 

(2.21) -W xx - d q W - U 2 W + dV + K = 

7 

and the boundary condition 

(2.22) W x (0) = W x (l) = 0. 

Let us recall that, for every 7 > 0, system (2.19)-(2.22) has the trivial solution 
(u ,v ,w ) = (0, 0,K O /d g ). Moreover, if Kq > 6, see (2.11), we have two other constant 
stationary solutions (u±,v±,w±) defined in (2.12). 

In the nonhomogeneous case, first, we consider U(x) and V(x), which are positive for 
all x e [0, 1]. Hence, under the assumption a > d c , using equations (2.19)-(2.20) we obtain 
the following identities 

(2.23) U(x) = ^V(x) and V(x) - ^ + d) 1 



d c y ' (a-d c y W(xY 

Moreover, adding equations (2.20) and (2.21), and using relations (2.23) we obtain the 
following boundary-value problem for the function W: 

(2.24) -W" - d g W - d b d } {dh+ J] -L + k = 0, 

7 (a-d c ) 2 W 

(2.25) W'{0) = W'{1) = 0. 
Here, to simplify our notation, we introduce the function 

(2.26) h(w) = -d 9 w - d fi {dh + d) 2 - + 



d c ) 2 w 



and we define number 



7T 2 



(2-27) 70 = T77— > 0, 

where W- is the constant from (2.12). 

First, let us summarize our results (proved in Section 5) on the nonexistence of solutions 
either to system (2.19)-(2.22) or to boundary value problem (2.24)-(2.25). 

• Under the assumption a < d c , the trivial solution (0, 0, K /d g ) is the only nonneg- 
ative solution of problem (2.19)-(2.22), see Proposition 5.1. 

• If a > d c and Kq < 0, then the boundary value problem (2.24)-(2.25) has no 
positive solutions, see Proposition 5.2.L Hence again, the trivial solution is the 
only nonnegative solution of (2.19)-(2.22). 

• If a > d c and = G, then the constant function w = Ho/d c is the only solution 
of (2.24)-(2.25), see Proposition 5.2.ii. 

• Let a > d c and > 0. The parameter 70 is the critical value of the diffusion 
coefficient in the sense that for 7 < 70, the boundary value problem (2.24)-(2.25) 
has no non-constant solutions., see Theorem 5.6. 
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Figure 2.1. All positive solutions of problem (2.24)-(2.25) with 7 > 
satisfying 3 2 7 < 7 < 4 2 7 . 

Hence, continuous patterns, i.e., non-constant continuous positive solutions, of system 
(2.19)-(2.22) can exist only if a > d c , Kq > 0, and 7 > 70. Here, we show that this is 
indeed the case. 

Definition 2.5. Let k E N and k > 2. We call a function W E C([0, 1]) a periodic 
function on [0, 1] with k modes ifW = W(x) is monotone on [0, ^] and if 



(2.28) 



W(x) 



W(x- 2 i) for 



x 



2j+2 



x) for x E [ 



c L fc ' fc J 
2.J+1 2j+2" 



for every j E {0, 1, 2, 3, ...} such that 2j + 2 < A;. 

We are now in a position to describe all possible nonnegative solution of the boundary 
value problem (2.24)-(2.25). 

Theorem 2.6 (Existence of continuous patterns). Assume that a > d c , Kq > 0, and 

7 > 70, where 70 is defined in (2.27). Consider the biggest n E N snc/i i/iat 7 > n 2 7 . 
T/ien, i/ie boundary value problem (2.24) -(2.25) has the following solutions: 

• the constant steady states w±, 

• a unique strictly increasing solution and a unique strictly decreasing solution, 

• for each k E {2, ...,n}, a unique periodic solution Wk with k modes that is increas- 
ing on [0, -r] as well as its symmetric counterpart: Wk(x) = Wk(l — x). 

There are no other positive solutions of the boundary value problem (2.24)-(2.25). 

All positive solutions of problem (2.24)-(2.25) with with 7 satisfying assumptions of 
Theorem 2.6 for n = 3 are sketched in Fig. 2.1. 

Remark 2.7. Using relations (2.23), we immediately obtain all solutions (U(x), V(x), W(x)) 
of system (2.19)-(2.22) such that U(x) > and V(x) > for all x E [0, 1]. 

Remark 2.8. It follows from the proof of Theorem 2.6 that non-constant solutions of 
problem (2.24)-(2.25) oscillate around the constant solution w_ and they are all below 
the constant w + . 
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In our construction of non-constant solutions of the boundary value problem (2.24)- 
(2.25), we use the well-known method from the classical mechanics that appears in the 
study of conservative systems with one degree of freedom (see e.g. the Arnold book [1]). 
In Section 5, we recall that approach for reader's convenience. Our main contribution to 
this theory consists in the fact that we have found the optimal value of the coefficient 7 
for which we have non-constant stationary solutions. Moreover, we prove that there is 
no other solutions, which is a consequence of the monotonicity of the function T = T(E) 
defined in (5.8), below. 

Our next goal is to discuss nonnegative solutions of (2.19)-(2.22) such that U(x) = 
V(x) = on a certain set I C [0, 1] which we refer to as a null set of a solution (U, V, W). 
Here, we do not expect W to be a C 2 -function and, in the following, we say that a vector 
(U, V, W) G L°°([0, 1]) x L°°([0, 1]) x C^QO, 1]) is a weak solution of system (2.19)-(2.22) 
if the algebraic equations (2.19)-(2.20) are satisfied for almost all x G [0, 1] and if 

(2.29) -- I W'{x)<f'{x)dx+ [ (-d g W(x)-U 2 (x)W(x) + dV(x) + K ) v (x)dx = 
7 Jo Jo 

for all if G ^([0,1]). 

Comparing to Theorem 2.6, the set of weak solutions is more complicated. 

Theorem 2.9 (Existence of discontinuous patterns). Assume that a > d c and Kq > O. 

There exists a continuum of weak solutions of system (2.19)-(2.22) with some 7 > 0. Each 
such a solution (U, V, W) G L°°(0, 1) x L°°(0, 1) x C 1 ([0, 1]) has the following property: 
there exists a sequence = x < X\ < x 2 < ... < x^ = 1 such that for each k G {0, N — 1} 
either 

• for all x G (xk,Xk+i), U(x) = V(x) = and W(x) satisfies -W" — d g W + k = 0, 



or 



• for all x G (xk,Xk+i), U(x) > and V(x) > are given by relations (2.23) , where 
the function W is a solution of equation (2.24). 

In Theorem 2.9, we do not attempt to classify all discontinuous stationary solutions (as 
it was done in the continuous case in Theorem 2.6), because they appear to be unstable 
solutions of the reaction diffusion equations (2.1)-(2.4), see the next subsection. Instead, 
in the proof of Theorem 2.9, we present a simple geometric argument which allows us to 
construct such discontinuous patterns. Analogous constructions, using either geometrical 
or analytical arguments can be found in [2, 9]. 

2.5. Instability of patterns. Finally, we discuss stability (in the sense of Lyapunov) of 
stationary solutions of system (2.1)-(2.4), which are constructed in Theorems 2.6 and 2.9, 
and surprisingly, we prove that they are all unstable. Let us be more precise. 

First, we consider a stationary solution (U(x), V(x), W(x)) of (2.1)- (2. 4) where W(x) is 
one of the functions from Theorem 2.6 and U(x) > 0, V(x) > are obtained from W(x) 
via relations (2.23). As the usual practice, we linearize system (2.1)-(2.3) around the 
steady state (U(x), V(x), W(x)) to obtain a system with three linear evolution equations 
with the linear operator C defined formally by 

(2.30) 
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where 

/ 4 (f - 1) ^1! o 

(2.31) A(x) = (a ij {x)) id=1 , 2 , 3 = 2K -d h -d ^ 

\ ~ 2K d ~ d 9 - WW) ) 

with the constant K = U(x)W(x) = dc ^^ , see the beginning of Section 6 for more 
details. We consider £ as a linear operator in the Hilbert space % = L 2 (0, 1) x L 2 (0, 1) x 
L 2 (0, 1) with the domain D{£) = L 2 (0, 1) x L 2 (0, 1) x W 2 ' 2 (0, 1). 

It follows from direct calculations that the constant coefficient matrix 

(2 ' 32 > A ^( dc{ l K 1] 

obtained from A(x) after removing the last row and the last column, has two real eigen- 
values of opposite sign. The positive eigenvalue Ao of A12 (given explicitly by formula 
(6.20)) plays the crucial role in the proof of instability of stationary solution to problem 
(2.1)-(2.4). 

Theorem 2.10 (Instability of continuous patterns). Consider the linear operator £ de- 
fined in (2.30), where W(x) be one of the functions from Theorem 2.6 and U(x), V(x) 
are obtained from W(x) via relations (2.23), except the constant solution (u + , v + , w + ), see 
(2.12). Then, the positive eigenvalue X of the matrix An belongs to the continuous spec- 
trum of the operator £. Moreover, there exists a sequence {A n } nS N of positive eigenvalues 
of the operator £ that satisfy X n — > X as n — > 00. 

Remark 2.11. For simplicity of the exposition, the constant stationary solution (tt+, v + , w + ) 
is excluded from Theorem 2.10, however, its instability is clear. Namely, it is an unstable 
solution of the kinetic system (A.1)-(A.3), see Corollary A. 10. 

Notice that, in order to present our idea in the simplest context, first, we prove Theo- 
rem 2.10 in the case of the constant stationary solution (U(x), V(x), W(x)) = (u_, v_, 
see Theorem 6.4, below. 



Corollary 2.12. Every continuous stationary solution (U, V, W) considered in Theorem 
2.10 is an unstable solution of the nonlinear system (2.1)-(2.4). 

This corollary results immediately from the classical theory. Indeed, by Theorem 2.10, 
every continuous stationary solution of system (2.1)- (2.4) is linearly unstable because the 
corresponding operator £ has infinitely many positive eigenvalues. Now, it suffices to 
note that — £ is a sectorial operator, see the monograph by Henry [3], hence, applying 
the general result from [3, Thm. 5.1.3] we show that (U, V, W) is an unstable solution of 
(2.1)-(2.4). 

Next, we show the instability of discontinuous stationary solutions from Theorem 2.9. 
First, we have to emphasize that a direct linearization of system (2.1)-(2.4) is not possible 
at a weak stationary solution satisfying U(x) = V(x) = at some x G [0, 1]. This is a 
consequence of the fact that the gradient of the nonlinear term in equation (2.1) is not 
well-defined when u = v = 0, see the Jacobian matrix (6.2), below. Hence, we modify our 
proof of instability in the following way. 

Let {Ux{x),Vx{x),W x {x)) be weak solution of (2.19)-(2.22) with a null set X C [0,1], 
namely, we assume that 

U x (x) = V x {x) = for x E X and U x {x) > 0, V x {x) > for x E [0, 1] \X. 
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Let us exclude the case X = [0, 1], where the constant vector (Uz, Vz, Wz) = (0, 0, no/d g ) 
is asymptotically stable solution of (2.1)- (2.4) by Corollary 2.4. For a null set X, we define 
the associate L 2 -space 

L|(0, 1) = {v G L 2 (0, 1) : v(x) = for x E X}, 

supplemented with the usual L 2 -scalar product, which is a Hilbert space as the closed 
subspace of L 2 (0, 1). Obviously, when the measure of X equals zero, we have L x (0, 1) = 
L 2 (0, 1), thus, we assume in the following that X has a positive Lebesgue measure and is 
different from the whole interval. 

Now, observe that if uq(x) = vq(x) = for some x G [0, 1] then by equations (2.1)-(2.2), 
we have u(x,t) = v(x,t) = for all t > 0. Hence, in view of Theorem 3.1, the space 

(2.33) Hz = L x (0, 1) x L|(0, 1) x L 2 (0, 1) 

is invariant for the flow generated by system (2.1)-(2.4) (notice that there is no "X" in 
the last coordinate of "Hi). The crucial part of our analysis is based on the fact that, as 
long as we work in the space Hz, we can linearize system (2.1)-(2.4) at a discontinuous 
steady state (Uz, Vz, Wz). Moreover, for each x G [0, 1] \X, the corresponding linearized 
operator agrees with £ defined in (2.30) with the matrix A(x) from (2.31). Hence, the 
analysis from the proof of Theorem 2.10 can be directly adapted to discontinuous steady 
states and we obtain the following counterpart of Corollary 2.12. 

Corollary 2.13. Every discontinuous weak stationary solution (Uz,Vz,Wz) of problem 
(2.19)-(2.22) with a null set X C [0,1], which was constructed in Theorem 2.9, is an 
unstable solution of the nonlinear system (2.1)-(2.4), considered in the Hilbert space Hz- 

Remark 2.14. In other words, the instability described by Corollary 2.13 appears when 
we perturb a stationary solution (Uz, Vz, Wz) on the set [0, 1] \X, namely, in those points 
where Uz and Vz are nonzero. On the other hand, by Theorem 2.2, small perturbations of 
this steady state on the null set X, leads to a solution satisfying u(x, t) — > and v(x, t) — > 
as t — y oo, for each x G X. 

The proofs of Theorem 2.10 and Corollary 2.13 are given in Section 6. 

2.6. Pattern formation. A natural question arises if our model exhibits a formation of 
any pattern, which persist for long times. Numerical simulations indicating the growth of 
such, computationally stable, spatially heterogeneous solutions were performed in [5, 6, 7]. 
The observed patterns take form of either periodic or irregular spikes and, now, it is 
clear that they cannot be described by stationary solutions, neither by smooth one from 
Theorem 2.6 nor discontinuous one provided by Theorem 2.9. Our present research is 
focused on understating these phenomena and answering questions on pattern formation 
in these kind of models. 

3. Existence of solutions 

We begin our study of properties of solutions to the initial value problem (2.1)-(2.5) by 
showing that it has a unique and global-in-time solution for all bounded and nonnegative 
initial conditions. 

Theorem 3.1. Assume that initial conditions Uo,v G L°°(0, 1) and Wq G jy 1,2 (0, 1) are 
nonnegative. Then, for every T > 0, the initial value problem (2.1) -(2.5) has a unique, 
global-in-time, nonnegative solution u G C([0, T], L°°(0, 1)), v G C([0, T], L°°(0, 1)) and 
w G C^^W 1 ' 2 ^,!)) such that u(x, •), v(x, •) G C\[0,T}) for every x G [0,1] and 
w G C\[0, T\- W^ 2 {0, 1)) n C((0, T}; W 2 > 2 {0, 1)). 
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Remark 3.2. We can improve regularity of solutions using results from [13, p. 112] in 
the following way. Assume that there exists a G (0,1) such that Uq,v G C q (0,1) for 
some a G (0,1) and if w G l7 2+q (0, 1) satisfies (w ) x (0) = (w ) x (l) = 0. Suppose, 
moreover, that uo(x) > 0, Vo(x) > 0, and wo(x) > for all x G [0, 1]. Then, the 
initial value problem (2.1)-(2.5) has a unique, global-in-time, and nonnegative solution 
u,v G C a < 1+a / 2 {[0, 1] x [0,T]) and w G C 2+a ' 1+a / 2 ([0, 1] x [0,T]) for every T > 0. 

The proof of Theorem 3.1 is more-or-less standard and we recall it for the completeness 
of the exposition. Below, we have only sketched the reasoning and we refer the reader to 
[13, pp. 108-123] for other details. 

First, we state a result on the local-in-time existence of solutions to the following initial- 
boundary value problem 

(3.1) in = fi(u,v) = ( -d c )u, 

\\u\ + \v\ J 

(3.2) v t = f2(u, v, w) = —dbV + u 2 w — dv, 

(3.3) w t = -w xx - d g w + f 3 (u, v, w) = -w xx - dgiu - u 2 w + dv + « , 

7 7 

supplemented with the homogeneous Neumann boundary conditions 

(3.4) w x (0, t) = w x (l, t) = for all t > 
and with initial conditions 

(3.5) u(x, 0) = u (x), v(x, 0) = v (x), w(x, 0) = w (x). 

Lemma 3.3. Let Uq,v G L°°(0, 1) and w G W 1,2 (0, 1). There exists T > such that 
the initial value problem (3.1)-(3.5) has a unique local-in-time solution 

u G C([0, T], L°°(0, 1)), v G C([0, T], L°°(0, 1)) and w G C([0, T], VT 1,2 (0, 1)). 

Proof. Here, it suffices to apply the abstract results by Rothe [13, Thm. 1, p. 111]. We 
recall that a mild solution (u, v, w) of problem (3.1)-(3.5), on a time interval [0, T) and with 
initial data uo,vo,wo G L°°(0,1), are measurable functions u,v,w : (0,1) x (0, T) — » R 
satisfying the following system of integral equations 

(3.6) t) — uq(x) + / fi(u(x, s),v(x, s)) ds, 

Jo 



(3.7) = v (x) + / f 2 (u(x,s),v(x,s),w(x,s))ds, 

Jo 

(3.8) w(x, t) = 5 , (t)u'o( a:; ) + / S(t — s)(f 3 (u,v,w))(x, s) ds, 



where S(t) is a semigroup of linear operators associated with the equation Z t = , -f~ 1 Z xx — 
d g Z on the interval (0,1), with the Neumann boundary conditions. Observe that the 
function fi = fi(u,v) in equation (3.1) is Lipshitz continuous if we put /i(0,0) = 0. 
Hence, to construct local-in-time solutions of system (3.6)-(3.8), it suffices to apply the 
contraction mapping principle based on the Picard iterations. Details of such reasoning 
can be found in [13, Thm. 1, p. 111]. □ 



Next, we show that the solutions are nonnegative. 



REACTION-DIFFUSION MODEL OF EARLY CARCINOGENESIS 



13 



Lemma 3.4. Let the assumptions of Lemma 3.3 hold true. Assume that (u, v, w) is a so- 
lution of problem (3.1)-(3.5) on a certain interval [0, T]. Ifuo, v anduiQ are nonnegative, 
then u(x,t) , v(x,t) and w(x,t) are nonnegative for all x £ [0,1] and t £ [0,T]. 

Proof. Step 1. First, notice that the inequality w(x,t) > for all x £ [0, 1] and t £ [0,T] 
implies v(x,t) > for all x £ [0,1] and t £ [0, T]. Indeed, using equation (3.2) with 
nonnegative w we obtain the inequality v t (x, t) > —(db + d)v(x, t), which implies v(x, t) > 
e - {db +d)T VQ ^ for all x e [ Q) i] anG i t £ [0, T}. 

Step 2. Without loss of generality, we can assume that w (x) > 0. Indeed, it suffices 
to replace w by w + e with e > and to use the continuous dependence of solutions of 
problem (3.1)-(3.5) on inital conditions. 

For wq(x) > 0, by the continuity of w (see Lemma 3.3), we obtain w(x,t) > for all 
x £ [0, 1] and for sufficiently small t > 0. Suppose there are xq £ [0, 1] and t > such 
that w(x,t) > for x £ [0,1] and t £ [0,to) and such that w(xo,to) = 0. Note that 
sup 0<0<to || ?x(t) || oo < M for some constant M > by Lemma 3.3 and v(x,t) > for all 
x £ [0, 1] and t £ [0,to] by Step 1. Hence, using equation (3.3) we obtain 

w t > —w xx — d„w — M 2 w + k . 

7 

In this inequality, all derivatives make sense, because, by a standard regularity theory for 
parabolic equation, we obtain 

w £ C\[0, t ], W 1 ' 2 ^ 1)) n C((0, t ], W 2 > 2 (0, 1)). 

Since /to > 0, by a comparison principle, we obtain w(x,t) > for all x £ [0,1] and 
t £ [0,t ], which contradicts our hypothesis w(xo,t ) = 0. 

Step. 3. Let us finally find the lower bound for u. As long as v — v(x, t) is nonnegative, 
it follows from equation (3.1) that u t (x,t) > —d c u(x,t). Integrating this differential 
inequality we arrive at the estimate u(x, t) > e~ dct u (x) for all x £ [0, 1] and t £ [0, t{). □ 

Proof of Theorem 3.1. Obviously, as long as a solution (u, v,w) of (3.1)-(3.5) is nonneg- 
ative, it is also a solution of our original problem (2.1)-(2.5). 

As the usual practice, to prove that local-in-time solutions from Lemma 3.3 exists for 
allt > 0, it suffices to show that for every T > 

SUp (|Kt)||oo + IK*) ||oo + !K*)l|oc) < OO, 
te[0,T) 

see eg. [13, Thm. l.iii, p. Ill] for more details. 

For positive u and v, we have v/(u + v ) < 1, hence, it follows from equation (2.1) 
that u t (x,t) < (a — d c )u(x, t) for all x £ [0, 1] and t £ [0, T). Integrating this differential 
inequality we obtain 

(3.9) sup ||u(t)||oo < e (a " dc)T ||n ||oo for every T > 0. 
te[o,T) 

Next, we use the integral equation (3.7) for v. Computing the L°°-norm and using (3.9) 
we obtain 

(3.10) |K*)IU<lh>o||oo+ / {(d b + d)\\v(s)\\ 00 + C(T)\\w(s)\\ 00 )ds, 

Jo 

where, following (3.9), we denote C(T) = e 2 (°~ dc ) T ||uo||^ . 
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In a similar way, computing the L°°-norm of the integral equation (3.8) and using the 
well-known properties of the semigroup S(t) we obtain 

(3.11) IK*)||oc < ||«>o||oo + «o+ / (d||v(a)|| 0O + C7(r)||«;(a)|| 0O )ds. 

Jo 

Finally, adding inequalities (3.10) and (3.11) and using the Gronwall lemma, we complete 
the proof that sup t6 r 0T ) (||w(t)||oo + |l w (*)||oo) < 00 f° r every T > 0. □ 

Proof of Theorem 2.1. Integrating both sides of equations (2.2) and (2.3) with respect to 
x G [0,1], adding the resulting formulas, and putting \i = mm{d g ,db}, we obtain the 
differential inequality 



d_ f 1 

dt 



(y(x, t) + w(x, £)) dx < —ji / (y(x, t) + w(x, £)) dx + k . 



Therefore, 



(3.12) J (v(x,t)+w{x,t))dx < ij (v (x) + w {x)) dx j e'^ + ^(1- e" M *), 

which implies (2.14) and the bound by hq/^i in (2.13), because v,w are nonnegative. 

To show the second bound in (2.13), notice that, for nonnegative u, v, we have ujiu + 
v) < 1, hence, it follows from equation (2.1) that Ut < — d c u + av. Integrating this 
differential inequality, we obtain that the functions U(f) = J u(x,t) dx and V(t) = 
J* Q v(x, t) dx satisfy 

(3.13) U(t) < e~ dct U{0) + a [ e~ Mt - s) V{s) ds. 

Jo 

Since V is a bounded function and d c > 0, we immediately obtain 

nt/2 ft/2 

/ e- Mt ~ s) V(s) ds < e~ dct/2 / V(s) ds ^ as t -»■ oo. 
Jo Jo 

On the other hand, we have 

f e- d ^~ s) V{s) ds < \ sup V{s) ) -j- (1 - e~ dct/2 ) , 
Jt/2 \se[t/2,t] J d c 

where lim^oo sup sg ^ / 2t i V(s) < Kq/ ji by (2.14). Hence, computing the limit superior as 
t — > oo of both sides of inequality (3.13), we complete the proof of (2.13). 

Next, let B(x, y, t) be the fundamental solution of the equation Z t = •y~ 1 Z xx — d g Z for 
x G (0, 1) and t > 0, supplemented with the Neumann boundary condition. It is well 
known (see e.g. [13, pp. 19 & 25]) that B(x,y,t) is nonnegative, 

(3.14) I B(x,y,t)dy < e~ dat for all x G [0, 1], t > 0, 
Jo 

and there exists a constant C > such that 

(3.15) / S^y^^^dy^Cr^e-^Hulli 

Jo 
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for each v G -^ 1 (0, 1), all x G [0, 1], and t > 0. Hence, equation (2.3) can be written as the 
following integral equation 

(3.16) w (x,t) = [ B(x,y,t)w (y)dy 

t r l 







JO 



+ / / B(x,y,t- s)(-u 2 (y,s)w(y,s) + dv(y,s) + K )dyds 



which, in view of the positivity of B, w, v, via inequalities (3.14), implies 

(3.17) |Ht)||oo<IKIUe-^ + d f [ B(x,y,t - s)v(y,s) dyds + ^(1 - e"^). 

Jo Jo dg 

Moreover, it follows from (3.12) and (3.15) that the second term on the right-hand side 
of (3.17) can be estimated as 

d / B(x,y,t - s)v(y,s)dyds 



JO 



<dC [ (t-sy^e-^-^Hs^ds 
(3.18) Jo t 

<<2C(IK||i + |KI|i) / (t- s )-V2 e -d s (t- s ) e -M S ds 







rt 

,«0 



+ dC ^ / (i - a )-V2 e -^(*-) (i _ e -/-) ds . 
/Wo 

It is easy to prove that limt_ ) . 00 Jjf (t — s) _1 / 2 e -da (* -s )e -|W ^ s = q an( j 

limsup / (t - s )-V2 e -<i s (*- 5 ) _ e -^) rfs = d -i/2 / S -V2 e -- ds 

t-voo Jo ^0 

Hence, the limit relation in (2.15) is immediately obtained from (3.17) and (3.18). □ 

Remark 3.5. By an inspection of the proof of Theorem 2.1, we obtain the following global- 
in-time estimate of w — w(x, t) : 

(3.19) IK*)||oo < IkolU + Cd ( 2 + -M ( + ||wo||i + — 



dgj \ fl 

where C and ft are the constants from (3.18). 

4. Large time behavior of solutions 

In this section, we prove preliminary results on the large time behavior of solutions of 
problem (2.1)-(2.5). Let us first prove that, under the assumption a < d c , each positive 
solution of system (2.1)-(2.3) converges exponentially towards the trivial steady state 

(0,0, Ko/dg). 

Proposition 4.1. Let a < d c . Assume that (u,v,w) is a nonnegative global-in-time 
solution of problem (2.1)-(2.5) corresponding to a bounded initial condition. Then, there 
exist positive constants C\, CV, o\, cr 2 dependent on the parameters in system (2.1)-(2.5) 
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and C*2 dependent also on ||wo||oo an d \\ v o\\oo, such that 

(4.1) < u(x, t) < u (x)e- {dc ~ a) \ 

(4.2) < v(x, t) < v (x)e- idb+d)t + C lU 2 (x)te 



(4.3) 



w(;t) 

a. 



(T\t 

K 



< C 2 e 

L°° 



-<J2t 



for all x G [0, 1] and t > 0. 



Proof. Since u and v are nonnegative, we have v/(u + v) < 1. Hence, we easily deduce 
from equation (2.1) the following inequality Ut(x,t) < —(d c — a)u(x,t), which implies 
estimate (4.1). 

Since Theorem 2.1 yields K w = sup <>0 ||u>(t)||oo < oo, then using equation (2.2) and 
estimate (4.1), we obtain the differential inequality 

(4.4) v t (x, t) < -{d b + d)v(x, t) + K w ul{x)e- 2{dc - a)t , 

which implies (4.2). 

To show the exponential convergence of w = w(x,t) towards Ko/d g stated in (4.3), it 
suffices to use the integral representation (3.16) of w, inequalities (4.1) and (4.2), and to 
follow estimates from the proof of Corollary 2.4, below. 

□ 

Proposition 4.2. Assume that (u(x,t),v(x,t),w(x,t)) is a nonnegative solution of prob- 
lem (2.1) -(2.5). Fix x G [0, 1]. The following two conditions are equivalent 

(i) u(x,t) — j- as t — )• oo and (ii) v(x, t) — >0 as t — > oo. 

Proof. We rewrite equation (2.2) in the following integral form 

v(x, t) = e-( db+d)t v (x) + / e-^+ d ^- s V(x, s)w(x, s) ds. 



By Theorem 2.1, the function w(x, t) is bounded for x G [0, 1] and t > 0, hence, combining 
condition (i) with the Lebesgue dominated convergence theorem in this integral equation, 
we can prove (ii). 

Recall that for nonnegative u and v we have u/{u + v) < 1. Hence, it follows from 
equation (2.1) that u t (x,t) < —d c u(x,t) +av(x,t). Integrating this differential inequality 
and following the previous argument we prove directly that (ii) implies (i). □ 

Before proving Theorem 2.2, in the following lemma, we show a kind of pointwise 
stability of the trivial steady state. 



Lemma 4.3. Let x G [0, 1]. Under the assumptions of Theorem 2.2, it holds 

(4.5) <u(x,t) < M + and 0<v(x,t) < M 

for all t > 0. 

Proof. Recall that the nonnegativity of u and v has been shown already in Lemma 3.4. 

Step 1. First, we suppose that there is T > such that u(x,t) < M(l + d c /a) for all 
t G [0, T}. Remember that < w(x, t) < K w for all x G [0, 1] and t > 0. Hence, it follows 
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from equation (2.2) satisfied by v that 

^ (v{x, t) - (^j A?j = - (4 + d) (v{x, t) - (^] M ) + u 2 w 



_(d 6 + d )^_£j M 

(4.6) 

< - (d b + d) [ v(x,t) _ ( ~ ) M 

M 2 (1 + ^V K w -(d b + d) (^Ym 



+ 

The term in the brackets on the right hand side of (4.6) is non-positive, because of 
assumption (2.17), hence, we skip it. Integrating the resulting differential inequality, we 
obtain 

(4.7) v{x, t) - (^j 2 M < e-^+W (v (x) - (^] ' M 

for all t G [0, T\. Using the assumption of vq(x) stated in (2.18), we see that the right-hand 
side of (4.7) is negative, hence, we obtain 

v(x,t)<(^j M for all t e [0,T]. 

Step 2. Now, suppose that there is T > such that v(x, t) < (d c /a) 2 M for all t G [0, T]. 
Applying the inequality u/(u + v) < 1 in equation (2.1) yields 

/ d \ 2 

u t (x, t) < av(x, t) — d c u(x, t) < a I — J M — d c ii(x, t). 

Integrating this differential inequality and using the assumption on uq(x) from (2.18), we 
obtain 

(4.8) u(x, t) < e-^uoix) + j (^j M (l - e"^) < M (l + ^ 
for all t G [0,T]. 

Conclusion. By Theorem 3.1, for each x G [0, 1], the functions u(x,t) and v(x,t) are 
continuous with respect to t. Suppose that inequalities (4.5) hold true for all t G [0,T) 
with some T > and at least one of them becomes an equality for t = T. Such a 
hypothesis, however, contradicts either the implication from Step 1 or from Step 2. This 
completes the proof of Lemma 4.3. □ 

Proof of Theorem 2.2. We are going to construct a sequence {T 2n } n >i, independent of 
x G [0,1], satisfying T 2n < T 2 ( n +i) an d \im n ^ 00 T 2n = +oo, such that u(x,t) and v(x, t) 
fulfill for each n G N the following inequalities 



(4.9) u(x, t) < 9 71 - 1 — 1 + — M and v(x, t) < 6 n — M for all t > T 2n . 
ci \ a i \ / 

Here 9 = [1 + (d c /a) 2 ]/2 < 1, because of the assumption d c < a. Since u and v are 
nonnegative, inequalities (4.9) imply lim^oo u(x, t) = lim^oo v(x, t) = 0. 
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We proceed by induction with respect to n. 

Step 1. By Lemma 4.3, we have v(x,t) < (d c /a) 2 M for all t > 0. The first inequality 
in (4.8) and the assumption on uq(x) lead to the estimates 

(4.10) u(x,t) < e-^u (x) + (-) M(l-e-^)<e- 4t M+-M(l-e- 4i ). 
Hence, choosing T\ > such that e~ dct M < (^) 2 M for all t > T x we obtain 

(4.11) u(x, t) < —M (1 + — J for all t > T x . 

a \ a J 

Consequently, it follows from the differential equation (2.2) and from the assumption on 
w(x, t) that 

2 



(4.12) Vt (x,t) < -(d b + d)v + 



d c , r ( \ d c 
—M 1 + — 



K w for all t > Ti. 



Integrating this differential inequality for t > T 1; we obtain 

(4.13) v{x,t) < e-^+W-^v&T!) + [a) } , (l - e'^^) 

db -\- d 

Using the assumption (2.17), written in the form 

{^) 2 M 2 (l + t)^ < (dX M 



db + d 

and the inequality v(x,T\) < (d c /a) 2 M from Lemma 4.3, we obtain from (4.13) that 



(4.14) v(x,t) < e -(*+«0(*- T O M + (^j M for all t > T x . 

Now, we choose T 2 > T x such that 

e -(*-w)(t-2i)^y M < i ^^y M _^y M j forall t > T2 . 

Therefore, the right-hand side of (4.14) is estimated as 

1 / . ( d r \ 2 \ ( d c \ 2 ( d c \ 4 ( <i c x 2 



"^^2r-(7j M+ \-a) M =\-a) m f ° ra " 

where = [1 + (d c /a) 2 }/2 < 1. Consequently, we have found T 2 > and we have proved 
the inequalities 

u(x,t) < — M + — 1 and v(a;,t) < f— 1 M# forall t > T 2 , 
a \ a J \ a J 

which correspond to (4.9) when n = 1. 

Sfep Now, we assume that (4.9) holds true for n G N and we prove it for n + 1. 
Equation (2.1) together with the inductive hypothesis for v(x,t) lead to 

(4.15) u t < av(x, t) - d c u(x, t) < a9 n (—} M - d c u(x, t) for all t > T 2n . 
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Integrating the differential inequality (4.15) for t > T 2 k, we obtain 

u(x,t) <e- d ^- T ^u(^T 2n ) + >f^VM(l- e -Mt-T 2n )\ 

d c \a J 
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(4.16) 



<e -d c (t-T 2n )Qn-l ( ^_ I l + 



M (l 



M 

d c (t-T 2n )^ 



for all t > T 2 „. 



Now, we choose T 2n+ i > T 2n in such a way that that following inequality holds 

e -d c (t-T 2n ) 6 n-ldc ( 1+ dc\ M <( ^ 9 n\ " M for ^ t > y 

ay ay \ a J 

which is equaivalent to the following one 

e -d c (t-T 2n )^ + (d c / a )) < {d c /a)8 n+1 for all t > T 2n . 
Hence, for t > T 2n +i we obtain from (4.16) that 



(4.17) 



d, 



u(x, t) < [ —6 n ) M + 6 n —M = —M6 n ( 1 + — 



1+& 



Since (1 + {d c /a)9 n )/{\ + (d c /a)) < 1, we have 



(4.18) 



d r 



d 



u(x, t) < —M6 n 1 + — for all t > T 2n+l . 



Next, integrating equation (2.1) and using (4.18), we obtain that the function v(x,t) 
satisfies for t > T 2n+ i the following inequality 

v(x, t) <e- {db+d){t - T ^ +l) v(x, T 2n+l ) 

d b + d K ' ' 

It follows from the assumption (2.17) that 



fcMfl" (! + £)] K w ^(d L 
db + d 



MB 



2/i 



and from the inductive hypothesis (4.9) that v(x,T 2n+ i) < 6 n (d c /a) 2 M. Hence, choosing 
T 2 ( n+ i) > T 2n+ i in such a way that 



-(d b +d)(t-T 2n+1 )nn ( ^\ M < 1 

W " 2 

or equivalently that 



d, 



d, 



-(d b +d)(t-T 2n+1 ) 



we deduce from (4.19) that 



< 



— e n M - — e 2n M 



d, 



for all t > T 2(n+ i) 



(4.20) v(x,t)< 



1 fd 



2 V a 



9 n M 1 



for all t > T 2(n+1) , 



for all t > T 2i n+i). 
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Since [1 + (d c /a) 2 6 n ]/2 < 6, inequality (4.20) implies 

(4.21) v(x,t) < (^) M6 n+1 for all t>T 2{n+1) . 

Consequently, inequalities (4.18) and (4.21) hold true for all t > T 2 ( n +i) and the proof 
of the inductive step is complete. □ 

Proof of Corollary 2.4- It follows immediately from Theorem 2.2 that u(x,t) — > and 
v(x,t) — > as t — t- oo uniformly in x G [0, 1]. To show that \\w(t) — K>o/dg\\oo — > as 
t — > oo we use the fundamental solution B(x,y,t) of the equation Z t = / ~f~ 1 Z xx — d g Z 
for x G (0, 1) and t > 0, supplemented with the Neumann boundary condition. Noting 
(3.14) and the positivity of all functions in the integral representation of w from (3.16) 
we obtain the inequality (c.f. (3.17)) 

(4.22) \\ w (t) - Ko/dgWoo < \\w - KoM?||ooe~ dfl ' + d / / B(x,y,t - s)v(y, s) dy ds. 

Jo Jo 

\2 



Since v(x,t) < (d c /a) M for all x G [0, 1] and t > by Lemma 4.3, we see 

(' [ B(x,y,t-s)v(y,s)dyds< (^] M [ 
Jo Jo \ a J Jo 



t/2 ,i fd\ 2 r t/2 

e- ds[t - s) ds ^0 as t -> oo. 



On the other hand, 



/ / B(x,y,t- s)v(y, s)dyds < [ sup v(x,s))—(l 

J t/2 Jo \s£[t/2,t] J tig 



e 



-dgt/2 



) 



It follows from Theorem 2.2 that the right-hand side above tends to as t — > oo. Conse- 
quently, we conclude from (4.22) that \\w(t) — Ko/d g \\oo — > as t — > oo. □ 

5. Construction of patterns 

In this section, we discuss stationary solutions of system (2.1)- (2. 4) (or equivalently 
of system (2.19)-(2.22)). First, we limit ourselves to nonzero U(x) and V(x) to obtain 
relations (2.23) as well as the boundary value problem (2.24)-(2.25). 

Let us begin with an arbitrary C 1 -function h : (0, oo) — > R. By the change of variables 

(5.1) x i — y Tx, where T = ^/-y, 
one can transform the boundary value problem 

-w" + h(w) = 0, w'(0) =w'(l) = 0, with xG (0,1), 

7 

into the problem 

(5.2) w" + h(w) = 0, xg(0,T), 

(5.3) w'(0) = w'{T) = 0. 

Together with equation (5.2), we consider the corresponding system of the first order 
equations 

(5.4) w' = z, z' = —h(w). 

Hence, solutions of equation (5.2) satisfying the boundary conditions (5.3) correspond to 
trajectories of system (5.4) satisfying z(0) = z(T) for a certain T > 0. Here, we should 
recall that system (5.4) in autonomous in the following sense: if (w(x), z(x)) is a solution, 
then (w(x + Xq), z{x + Xq)) is a solution for all xq G M, as well. 
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H(w + ) 




E 



H = H(w) 



H(w-) - 



E W- W 2 ,E W + 



Z 




z = ±y/2{E - H(w)) 



W 



Figure 5.1. The graph of the potential energy H = H(w) defined in (5.6) 
and the closed trajectory of system (5.4) with the energy level E > 0. 



Multiplying second equation in (5.4) by z = z(x) and using first one, we obtain the 
equation for all trajectories of system (5.4) 



where H' = h and E e IR is an arbitrary constant. Recall that the constant E is called 
the total energy in the classical mechanics, the function H corresponds to the potential 
energy, and relation (5.5) is the first integral of system (5.4). 

It follows from equation (5.5) that all trajectories of system (5.4) are symmetric with 
respect to the w-axis. Hence, the condition z(0) = w'(0) = z(T) = w'(T) is satisfied for a 
certain T > if equation (5.5) describes a closed curve on the u>z-plane for a some £gR. 
Such closed curves exist only if the potential energy H = H(w) has a local minimum at 
a certain point see Fig. 5.1. 

Let us apply these classical ideas to obtain preliminary results on the nonexistence of 
stationary solutions of system (2.19)-(2.22). Now, the function h is defined in (2.26). 

Proposition 5.1. Assume that a < d c , 7 > 0, and d, db, d g , kq be arbitrary and positive. 
Then, the only nonnegative solution of system (2.19)-(2.22) is the trivial steady state 



Proof. If a < d c , the nonexistence of stationary solutions, except the trivial one, is guar- 
anteed by Proposition 4.1. 

For a = d c , it follows from equation (2.19) that U = and, by equation (2.20), we obtain 
V = 0. Consequently, using (2.21) we obtain the following boundary value problem for 



(5.5) 



2 



+ H(w) 



E. 



(0,0, Ko/dg). 
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the function W = W(x) 

-W xx -d g W + K = 0, x £ (0, 1) 

7 

W x (0) = W x (l) = 0. 

Here, the potential energy H(w) = —d g w 2 /2 + hqw has no local minimum and the cor- 
responding equation for trajectories (5.5) does not describe any closed curve, except 
E = H(K /d g ) and the stationary point (w, z) = (K /d g , 0). □ 

In the following proposition, we discuss the nonexistence of positive solutions of the 
boundary value problem (2.24)-(2.25). 

Proposition 5.2. Assume that a > d c and 7 > 0. Let other constants in (2.24) be 
nonnegative and arbitrary. 

i. If Kq < (cf. (2.11)), then the boundary value problem (2.24)-(2.25) has no 
positive solutions. 

ii. // Kq = 6 ; then the constant function w = Ko/d g is the only solution of the 
boundary value problem (2.24)-(2.25). 

Proof. Direct calculations (see Appendix A - a discussion around equation (A. 7)) show 
that, under the assumption Kq < O, the function h defined in (2.26) satisfies h(w) = 
H'(w) < for all w > 0. Hence, the associated potential energy H, (see formula (5.6), 
below) is strictly decreasing and, in consequence, equation (5.5) does not define any closed 
trajectory of system (5.4). 

If Kq = 0, the function H is decreasing, as well. Here, however, system (5.4) has a 
stationary point (n /(2d g ), 0). □ 

Hence, in view of Propositions 5.1 and 5.2, we have to assume that a > d c and k,q > G 
to be able to construct nonconstant stationary solutions of the boundary value problem 
(2.24)-(2.25). It follows from direct calculations (see the beginning of Appendix) that, 
under these assumptions, the numbers w± from (2.12) are the only positive zeros of 
the function h defined in (2.26). Moreover, the corresponding potential energy (namely, 
H' = h) 

d 2 c (d b + d) 
~dj 

has a local minimum in u>_, a local maximum in w + . 

To describe inhomogeneous stationary solutions, let us come back to the abstract 
boundary value problem (5.2)-(5.3) and recall how to calculate the number T > such 
that w'(0) = w'{T) = 0. Indeed, in view of equation (5.5), every solution w = w(x) to 
problem (5.2)-(5.3) satisfies the first order differential equation 

(5.7) w'{x) = ±y / 2(E - H(w(x))). 

Hence, choosing the upper branch in (5.7) and integrating with respect to x we obtain 

w'(x) dx 



(5.6) H(w) = -d g w 2 /2 - 4 7„ Tn 2 1o S w + ^ 



y/2(E - H(w(x)))' 

Next, let us introduce w^e < W2,e such that w^e = w(Q) and W2,e = u;(T), hence 
H(w\,e) = H{w2,e) = E, see Fig. 5.1. Since w = w(x) is nondecreasing for x £ [0, T], 
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by the change of variables y = w(x) we obtain the following formula for T as a function 
of E: 

dy 



(5.8) T = T(E) 

! W1<B ^2{E-H{y)) 

Our classification of all nonnegative solutions of the boundary value problem (5.2)- 
(5.3) is based on detailed analysis of the integral (5.8), which properties are stated in the 
following lemmas. 

Lemma 5.3. Assume H £ C 2 (0, oo) has a local minimum at W- such that H"{wJ) > 0. 
We consider every E > H{wJ) such that there exist and w 2j e with the following 
properties (see Fig. 5.1): 

• H(w ljE ) = H(w 2>E ) = E, H'(w 1>E ) ± 0, H'(w 2>E ) ^ 0, 

• W 1)E <W-< W 2 ,E, 

• H{w) < E for all w £ (wi jE , W2,e)- 

Then, for all such constants E the integralT(E) defined in (5.8) is convergent and depends 



continuously on E. Moreover, \im E \ iH ^_^T(E) = n/ ^H"(w-). 

Proof. By the assumptions on the function H, the integrand of T(E) from (5.8) has 
singularities at W\ >E and wi iE , only. Thus, the integral T(E) is convergent by the Taylor 
expansion, because H'(wi tE ) 7^ and H'{w 2 ,e) ^ 0. To show the continuous dependence 
of T(E) on E, it suffices to apply e.g. the Lebesgue dominated convergence theorem. 

To calculate the limit \im E ^H(w-)T(E), first, we recall that E = H(w 2jE ) and we 
consider the integral 

f W2 ' E dy r^E-w. d 

Ti{E) ~ 



^2(H(w 2 , E ) - H(y)) Jo y/2(H(w 2 , E )-H(y + w^y 

Defining the new parameter s = w 2 ^e — w~ and the shifted function H(w) = H(y + W-), 
moreover, changing variables we obtain 

(5.9) lim T X (E) = lim [' dy = lim / dy 

E-+HHD-) s^oJ Q ^J 2 (H(s)-H(y)) ^2s^ (H(s) -H(sy)) 

Consequently, we apply the Lebesgue dominated convergence theorem combined with the 
l'Hospital rule to the integral on the right-hand side of (5.9) to show 

lim Tl {E) = -=1= f -fi=. 

Jh'\o) Jo V 1 - y 2 

Now, it suffices to recall that h' (0) = H"{wJ). 

In a completely analogous way, one calculates the limit 

dy 1 f° dy 



lim 



E^Hfr-) J W1E y/2(E - H{y)) y/H"{w_) J-i s/Y^tf 

and the proof is complete because 1/ \/\ — y 2 dy = ir. □ 

Lemma 5.4. Assume that H £ C 2 (0,oo) has, for certain < u>_ < w + , the following 
properties 

• H is strictly decreasing on (0,iD_) and strictly increasing on (w_,w + ); 

• H'(w+) = 0. 
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For every E G (H(w_), H(w + )) , choose w^e < w~ < w 2< e < w + such that H(wi tE ) = 
H(w 2> e) = E, see Fig. 5.1. Then ]xm E *H(iv+)T(E) = +00. 

Proof. It suffices to show that 

f W2 ' E dy 
am I — = = +00. 

Jw. ^2{H{w 2 ,e) ~ H(y)) 

The function H is increasing on (w-, w+) and H'{w+) = 0, hence, by the Taylor expansion, 
we obtain 

777 

H(w 2tE )-H(y) <H(w + ) -H(y)<-(w + -y) 2 with m= sup \H"{y)\ 



for all y G (u>_ , w 2> e) ■ Consequently, 



dy 1 f W2 ' E dy 

> —= I — > +OO if W 2> E — > U)+ 



^2(H(w 2yE ) - H(yj) Vm w + - y 

□ 

Lemma 5.5. Assume that a function H has all properties stated in Lemma 5.4- Suppose, 
moreover, that H G C 3 (0,oo), h = H' satisfies h"(w) < for every w G (0,w + ), h" is 
not constant on any interval, and h(wJ) = H'iwJ) = 0. Then T{E) defined in (5.8) is a 
strictly increasing function of E G (H(w_), H(w+)). 

Proof. Obviously, it suffices to show that dT(E)/dE > using the explicit formula (5.8). 
These involved calculations were done by Loud [4] by using a clever change of variables. 
Let us recall that result in our particular case. 
First, we consider the integral 

" W2 > E dy 



(5.10) Ti(E) 

y/2(H(w 2 , E ) - H(y)) 

where shifting the variable y as in the proof of Lemma 5.3, we can assume that u)„ = 0, 
without loss of generality. Since H is strictly increasing on (0, w+), denoting H(A) = E, 
it suffices to show that the following function of A 

f A dv 

T(A) = Ti(H(A)) 



'2{H{A)-H{y)) 
is strictly increasing. It is calculated in [4, Thm. 1, Eq. (2.5)] that 



H{y)h\y) 



dy 



%) 2 2j y/2(H(A) - H(y))' 



Now, we denote by f(y) the quantity in the brackets on the right hand side of (5.11) and 
we observe that, for every y G [0, A], we have 

(5.12) ^-{h 2 {y)f{y))=H{y)h"{y) and Zi 2 (0)/(0) = 0. 

By the assumptions, H(y) > and h"(y) < 0. Therefore, it follows from (5.12) that 
h 2 (y)f(y) < and, in consequence, f(y) < for all y G [0,A]. Finally, the right-hand 
side of equation (5.11) is positive because h(A) > and H(A) > 0. 

The analysis is completely analogous in the case of the counterpart of the integral 
Ti(E), where we intergrate with respect to y G [w^BjW-]. □ 
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Now, we come back to boundary-value problem (2.24)-(2.25) with the function h defined 
in (2.26), and with the corresponding potential energy H from (5.6). 

Theorem 5.6. Assume that a > d c , k,q > G, and 7 G (0,7o] ; whereto is defined in (2.27). 
Then, the constant steady states w± (cf. (2.12)) are the only solutions of the boundary 
value problem (2.24)-(2.25). 

Proof. Recall that every non-constant solution w = w(x) of problem (2.24)-(2.25) corre- 
sponds to a non-constant trajectory (w(x), z(x)) of system (5.4) such that z(0) = z{T) = 
with T = J^f. However, by Lemmas 5.3 and 5.5, such trajectories exist only if T > 

n/y/H"(w..) = V?^. □ 

Now, we prove our main result on the existence of continuous and positive stationary 
solutions of system (2.1)- (2. 4). First, we notice that, if W — W(x) is a solution of 
problem (2.24)-(2.25), then so is W(x) = W(l — x). More generally, for every A,BeM, 
the function W(x) = W(Ax — B) is a solution of equation (2.24) with 7 replaced by A 2r y, 
satisfying the boundary conditions W X {B/A) = W X ((T + B)/A) = 0. 

Proof of Theorem 2.6. The potential energy H = H(w) associated with equation (2.24) 
has the form (5.6) and, under the assumption Kq > O, it has a local minimum in w_, 
a local maximum in w + , and all properties required in Lemmas 5.3-5.5. Obviously, the 
numbers w± are constant solutions of the boundary value problem (2.24)-(2.25). 

For every E G (if («)_), H(w + )), let us consider the trajectory (w(x), z(x)) of system 
(5.4) such that z > and z(0) = z(T) = (namely, the upper half of the trajectory drawn 
in Fig. 5.1). Here, T = T(E) is defined by the integral (5.8) and, by Lemmas 5.3-5.5, 
this is a continuous and increasing function of E, which takes all values from the half-line 
(n/ \J H"[wJ), 00). Due to the change of variables (5.1), the function W(x) = w(x/T) is 
the unique increasing solution of the boundary value problem (2.24)-(2.25) with 7 = T 2 . 
On the other hand, the lower half of the trajectory drawn in Fig. 5.1, more precisely, the 
function W(x) = W{1 — x), is the unique decreasing solution of (2.24)-(2.25). 

Now, for fixed 7 > 70, let us choose the biggest n G N such that 7/n 2 > 70 and suppose 
that n > 2. For each k G {2,...,n}, we consider the unique trajectory (w fc (x), Zk(x)) 
of system (5.4) such that z k (x) > and Zfc(0) = Zk(T k ) = with T k = y/j/k. Hence, 
Wk(x) = Wk(x/ (kTk)) for x G (0, kT k ) is the solution with k modes of the boundary value 
problem (2.24)-(2.25). Another solutions with k modes is the symmetric counterpart 
W k {x) = W(l-x). 

These are all solutions and there exists no other solutions, because T = T(E) is a 
continuous and strictly increasing function of E, see Lemmas 5.3-5.5. Hence, we obtain all 
trajectories (w(x), z(x)) of system (5.4) such that z(0) = z(T) for each T > n/ ^jH"(wJ). 

□ 

Proof of Theorem 2.9. First, we change variables as in (5.1), hence, we consider weak 
solutions of system (2.19)-(2.22) (with 1/7 = 1) on the interval [0, T] (with T = y/^y). In 
these new varables, if U(x) = V(x) = for some x G (0,T), the function W satisfies the 
equation 

W"(x) - d g W(x) + k = 0. 
Trajectories of the corresponding system on the iyZ-plane 



(5.13) 



W' = Z, Z' = dgW - K 
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Figure 5.2. Two closed trajectories of system (5.4) with energy levels 
E\ > 0. Three trajectories of system (5.13) with energy levels E 2 > 0. The 
bold line shows a trajectory corresponding to a discontinuous pattern. 

are unbounded and have the following explicit form Z = ±^2E 2 + d g W 2 — 2k W for 
every constant E 2 > 0, see Fig. 5.2. 

On the other hand, if U (x) ^ and V(x) ^ for some x £ (0, T), we use relations (2.23) 
for U(x), V(x) together with equation (2.24) for W(x) as well as the coresponding system 
on the VFZ-plane (5.4) with h defined in (2.26). Recall that, under our assumptions, this 
system has closed trajectories described by the formula (5.5), see Fig. 5.1. 

Now, we construct a continuous trajectory (W, Z), which corresponds to a weak solution 
of (2.19)-(2.22) in the following way. We begin at the IF-axis at x = and finishing at the 
W certain T > going along either trajectories of system (5.13) or trajectories 

of system (5.4). At each point of the intersection of two trajectories of different types, 
the function W is C l because Z = W is continuous. Fig. 5.2 shows examples of such 
trajectories. 

One can easily check that the constructed-in-this-way function W = W(x) satisfies the 
integral equation (2.29) for every test function tp £ C 1 ([0, 1]). □ 



6. Linearization and spectral analysis 

In this section, we show instability of stationary solutions of system (2.1)-(2.5), which 
are constructed either in Theorems 2.6 or in 2.9. Writing equations (2.1)-(2.3) in the form 

(6.1) u t = fi(u,v), v t = f 2 (u,v,w), w t = -w xx + f 3 (u,v,w), 

7 

we obtain that the differential of the mapping F = (fi, f 2 , fs) : K 3 — > M 3 satisfies 

-d r . ^ 



(u+v) 2 c (u+v) 2 

(6.2) DF(u, v, w) — | 2uw -d b -d u 2 

-2uw d —d g — u 2 

Assume first that U(x) > and V(x) > for some x £ (0, 1). Using relations (2.23), 
we obtain the equalities 

V(x) 2 _ d\ ^ U{xf _{a- d c f 



(U(x) + V(x)) o 2 (u( kX ) + V { x )) o? 
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Hence, 

/ dc(a-dc) (a-dc) 2 q 

(6.3) DF(U, V, W) = A{x) = (a,,),., ,. 2 . :i = 



K 2 

w 2 (x) 



2K -d b - d 



with the constant 

v TU Sixr , . d c (d b + d) 

K = U(x)W(x) = — — L , 

a — d c 

see (2.23). Notice that only the coefficients a 23 and 033 depend on x. 

Hence, the linearization of problem (2.1)-(2.5) at a continuous steady state (U,V,W) 
(namely, when U(x) > and V(x) > for all x G [0, 1]) contains the linear operator 

(<p\ /o 

(6.4) £^ = 00 

V v J V ±d 2 v 

and we consider it as an operator in the Hilbert space H with the domain D(C), where 

(6.5) U = L 2 {0, 1) x L 2 (0, 1) x L 2 (0, 1) and D(C) = L 2 {0, 1) x L 2 (0, 1) x W 2 ' 2 {0, 1). 

First, we prove two lemmas which allow us to characterize the point spectrum of the 
operator C Here, we do not use any particular form of the coefficients in the matrix A 
in (6.3). Below, we denote by A12 the matrix obtained from A after removing the third 
row and the third column, namely, 




(6.6) A 



12 



an 012 

Q<21 a 22 



Lemma 6.1. Assume that the matrix A12 has constant coefficients and let Ao be its 
eigenvalue. Then A belongs to the continuous spectrum of the operator C defined in 
(6.4). 

Proof. Let us first show that Ao is not an eigenvalue of the operator C To do it, we should 
show that (<f, if}, 77) = (0, 0, 0) is the only solution of the system 

(an - A )y2 + a i2 ^ = 0, 

a 2 if + (a 22 - A )^ + a 23 ?7 = 0, 



(6.7) 



-77" + a 3 iv9 + a 32 ip + (a 33 - A )?7 = 0, 



7 



r/(0) = t}'(1) = 0. 

Since A is an eigenvalue of the matrix A12, the vectors (a u — A , a i2 ) and (a 2 i, a 22 — A ) 
are linearly dependent, namely, there is r\ G R such that 

(6.8) (an - A , a 12 ) = r 1 (a 2 i, a 22 - A ). 

Hence, it follows from the first and the second equation in (6.7) that rj = and there 
exists a number r (in fact, r = — ai 2 /(a n — A )) such that (p = rip. We substitute this 
relation into the third equation in (6.7) to obtain the equality (ra3i + 032)^ = that 
implies ip = and, consequently, <p = 0. 

Thus, we have proved that the operator C — \qI : D(C) — > H (cf. (6.5)) is invertible. 
However, its inverse cannot be continuous, because (C — XqI)^ 1 is not defined on the 
whole space "H. Indeed, using relation (6.8) in the first and the second equation in (6.7), 
we obtain the equality rxa^r/ = rig — f. Hence, 

(6.9) (r ig -f)/( ri a2 3 )eW 2 ' 2 (0,l) 
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is a necessary condition for system (6.7) to have a solution (ip,ip,r]) G D(C). Obviously, 
the condition in (6.9) is not satisfied for every f,gE L 2 (0, 1). This completes the proof 
of Lemma 6.1. □ 

Lemma 6.2. A complex number X is an eigenvalue of the operator C if and only if the 
following two conditions are satisfied 

• X is not an eigenvalue of the matrix A\ 2 , 

• the boundary value problem 

1 .. detL4-AJ) , . 

-V + i— r4 7T\V = 0, x G 0, 1 

(6.10) 7 det(A 12 -XI) 1 v ; 

r/(0) = rj'(l) = 
has a nontrivial solution. 

Proof. Assume that the number A G C is an eigenvalue of the operator C. Hence, the 
system 

(an - A)y? + a 12 i/j = 

(6.11) a 2 i(f + (a 22 - X)ip + a 23 r] = 
\rf' + a 3lV? + a 32^ + (ct33 - A)?7 = 0, 

supplemented with the boundary condition 7/(0) = ^'(1) = 0, has a non-zero solution 
(ipo, "00) ^o)- By Lemma 6.1, A is not an eigenvalue of the matrix Ai 2 . Hence, det(^4i2 — 
XI) and from the first and the second equation in (6.11) we obtain 

/^ 10 n Oi 2 a 2 3 j . («n — A)a 2 3 

(6-12) if = - — i—. —prt] and w = — - — ttt 7 ?- 

v ' v det(A 2 -XI)' v det(^li 2 - XI) ' 

Substituting these identities to the third equation in (6.11) provides the relation 

/ C1 „, l „ , 031012023 - 032(011 - A)a 2 3 + (033 - A) det(.4i2 - XI) 

6.13 -77 + — — — r/ = 0, 

7 det(^4i 2 - XI) 

which reduces to the equation in (6.10), because the quantity 031012023 — 032(011 — A)a23 + 
(a 33 — A) det(^4i2 — XI) is the Laplace expansion of the determinant det(^4 — XI) with 
respect to its third column. 

Now, let us prove the reverse implication. Assume that A G C is not an eigenvalue of 
A± 2 and denote by 77 a non-zero solution of problem (6.10). Hence, the vector (ip , ip , 77) , 
where ip and ip are defined in (6.12), is the eigenvector of the operator £ corresponding 
to the eigenvalue A. □ 

We are now in a position to show the existence of an infinite sequence of positive 
eigenvalues of the operator L and to prove Theorem 2.10. First, we consider the simplest 
case of the constant stationary solution (U (x), V(x), W(x)) = (u_, u_, tD_). Here, we need 
a simple technical lemma. 

Lemma 6.3. Assume that a 3 x 3-matrix A with real coefficients has all 3 eigenvalues 
with negative real parts. Then, its characteristic polynomial satisfies det(*A — XI) < for 
all A > 0. 

Proof. Since the polynomial det(A— XI) has real coefficients, it has three roots Ai, A2, A3 G 
C satisfying Ai G R and A3 = A 2 . Hence, 

det(A-XI) = -(A-Ai)(A- A 2 )(A-A 2 ) = -(A - Ai)(|A 2 | 2 - 2ARe A 2 + A 2 ). 
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The factor (A — X\) is positive for all A > 0, because Ai < by the assumption. The last 
factor on the right-hand side is positive for every A > 0, because ReA2 < 0. □ 

Theorem 6.4 (Instability of the constant steady state). Denote by Ao the positive eigen- 
value of the matrix A\%. Consider the operator £ from (6.4) with the constant coefficient 
matrix A obtained from the constant steady state (w_, w_, w_). Then, there exists a se- 
quence {\ n }n£N of positive eigenvalues of the operator C that satisfy A n — > Ao as n — > oo. 

Proof. In view of Lemma 6.2, we look for an infinite sequence of numbers X n such the 
boundary value problem (6.10) has a nontrivial solution. First, recall that the following 
eigenvalue value problem on the interval [0, 1] 

(6.14) -r]" + nr) = 0, t/(0) = t/(1) = 

7 

has a nonzero solution rj(x) = r] n (x) = cos(nirx) if fi = fi n = n 2 7r 2 /7 for each n G N. 

By Lemma 6.3, the polynomial det(A — XI) is negative for all A > 0, because the steady 
state (w_, w_) is a stable solution of the kinetic system (A.1)-(A.3), see Corollary A. 10 
in Appendix. 

Since Ao is the only positive eigenvalue of the matrix A±2, the quadratic polynomial 
det(Ai2 — XI) has the following properties: det(A±2 — XI) < for A G [0, Ao) and det(Ai2 — 
X I) = 0. Hence, det(^4- A/)/det(^i 2 - XI) > for all A G [0, A ) and the left-hand side 
limit satisfies 

det(A-XI) 

lim — — — = +oo. 

a^a„ det(yli2 - XI) 

By continuity, we immediately obtain a sequence A n — > Ao such that 

det(.4-A re I) _ n 2 vr 2 



det( v 4i2 - X n T) 7 
for all sufficiently large n G N. □ 

More-or-less similar idea is used to find positive eigenvalues of the operator £ in the 
case of non- homogeneous steady states. First, however, we recall properties of eigenvalues 
\i G R of the boundary value problem 

(6.15) ~r)" + nqr) = 0, rf(0) = r/(l) = 0, 

7 

where q G L°°(0, 1) is nonnegative. It is well-known that problem (6.15) has a sequence 
of eigenvalues {/^(g)}^ satisfying 

= Ho(q) < /ii(g) < fi 2 {q) < fJa(q) < ► +oo. 

Moreover, these eigenvalues depend continuously on the potential q = q(x) in the following 
sense. 

Lemma 6.5. Let qi,q2 G L°°(0, 1) be nonnegative and fixed. Denote by {/-tn(<Zi)}£!Li, 
i = 1,2, the corresponding eigenvalues of problem (6.15). Then, there exists a constant 
C > such that for all n G N, it holds 

\Vn(q\) - A*n(S2)| < lC\\qi - g 2 ||L2( ,i). 

Moreover, if q\(x) < q 2 (x) a.e. then jJ> n (q 2 ) < H> n (Qi) f or eac h n G N. 
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Proof. Denoting by G(x, y) the Green function of the operator —d 2 /dx 2 with the Neumann 
boundary conditions, we convert problem (6.15) into the integral equation 

v( x ) = ta*(?) / G ( x ^y)i(y)v{y)dy, xe[o,i], 

Jo 

where the right-hand side defines a compact, self-adjoint, nonnegative operator on L 2 (0, 1). 
Hence, properties of eigenvalues, stated in this lemma, result immediately from the ab- 
stract theory, see eg. [12, Cor. 5.6 and Thm. 5.7]. □ 

Proof of Theorem 2.10. By Lemma 6.1, the number Ao belongs to the continuous spec- 
trum of the operator L. Next, denoting 

(6 ' 16) g(x ' A)= detCAu - A/) ' 

we rewrite problem (6.10) in the form 

(6.17) - rj"(x) + q(x, X)rj(x) = 0, r/(0) = r/(l) = 0. 

7 

By Lemma 6.2, it suffices to find a sequence A = A„ — > Ao such that problem (6.17) has a 
nonzero solution. We shall proceed in a sequence of steps. 

Step 1. There exists e > such that for every A G (Ao — s, Ao), we have 

min q(x, A) = lo(X) > and ^(A) —> oo and A — > A . 

xe[o,i] 

To show this property of the potential q(x, A), we use its expanded form from equation 
(6.13) 

(6.18) q(x, A) = Q2 3( x ) a31 ^ e 2 t( -^ 2 ^ Ql A 1 j) A ^ + Q 33( x ) ~ A ' 
where, by (6.3), 

(6.19) a 23 = \ and a 33 (x) = -d g - , with K = — , 

W 2 (x) W 2 (x) a - d c 

and other coefficients on the right-hand side are x- independent. Recall (see the proof of 
Theorem 6.4) that det(^4i 2 - XI) < for A G (0, A ) and det(A 12 - XqT) = 0. 

Next, we show that the numerator in the fraction on the right-hand side of (6.18) 
satisfies R(Xo) = 031012 — 032(0-11 — Ao) < 0. Indeed, using the explicit form of the 
coefficient of the matrix A12, we obtain 
(6.20) 

1 /_ ( dja - d c ) +db + d \ + l( d e (a -d c ) +db + d Y + 4 4(a - d c )(d b + d) 



2 \ \ a J V V a J a 

and 

^ (Aq) = _ 2 d c (a-d c )(d b + d) + d f d c (a - d c ) + Aq 
a \ a 

Hence, denoting y = (d c (a — d c ))/a and D = d b + d leads to 

d (y - D + y/(y + Dy + Wy) 
R(X ) = R(y) = -2Dy + '- 
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Obviously, R(0) = 0. We leave for the reader to check that dR/dy < for all y > 0, 
which implies that R(y) < for all y > 0. This implies immediately that R(\o) < 0. 
Hence, provided e > is sufficiently small, for all A £ (A — e, A ), we obtain 

i?(A) n , i?(A) . . 

> and - — -> oo as A -»■ A . 



det(.Ai2 - A/) det(A2 - A/) 

Since min xe [ 0i i] a 23 (x) > (see (6.19)) and since a 33 (x) is a bounded function, choosing 
smaller e > (if necessary) we obtain min a . g [o i i] <?(£, A) = a; (A) > for all A £ (A — e, A ) 
and u;(A) — > oo for A — > A . 

S'tep By the Sturm-Liouville theory and Step 1, for every A £ (A — s, A ), there 
exists an increasing sequence 

< m(q(-, A)) < fi 2 (q(-, A)) < ... < H n (q(-, A)) < ... +oo 
of eigenvalues of the problem 

-r/' + A)t? = 0, x£(0,l) 
(6.21) 7 

t/(0) = r/(l) = 0. 

Our goal is to show that there exists A„ £ (Ao — £, Ao) such that [i n (q(-, A n )) = 1. Then, 
the corresponding eigenfunction of (6.21) will be a non-zero solution of (6.17). 

Step 3. For each n £ N, the quantity /i n (g(-, A)) is a continuous function of A. Indeed, 
this is an immediate consequence of Lemma 6.5, because q(x, A) is a continuous function 
of A £ (A -£,A ). 

Step 4- Let us show that yU„(A) — > when A — > Ao- By Step 1, we have cu(A) = 
min^o,!] q(x, X) — > oo when A £ (Ao — s, Ao) and A — > Ao- Since u(X) is a positive 
constant for each A £ (Ao — s, Ao), the eigenvalues /i of the problem 

-rj" + fiu{X)r] = 0, t/(0) = r/(l) = 0. 

7 

are given explicitely /i n (u;(A)) = ri 2 7r 2 /(u;(A)7). Since w(A) < g(x, A), the comparison 
property for eigenvalues (stated in Lemma 6.5) implies 

Hn{q{-, A)) < fi n (u(X)) = — — ^0 as A -> A . 

w(Aj7 

Conclusion. Since < /ii(g(-,A)) < fi2(q(-,X)) < ... < /i n (g(-,A)) < ... — > +oo and 
A n (g(-, A)) — )• as A — )• A , for every n £ N, it follows from the continuous dependence of 
yU n (A) on A (see Step 3) that there exists A n — > A such that /i n (A n ) = 1, provided n is 
sufficiently large. This completes the proof of Theorem 2.10. □ 

Proof of Corollary 2.12. This is a consequence of Theorem 2.10 and of a general result 
from [3, Thm. 5.1.3] as it is explained in the paragraph following Corollary 2.12. □ 

Proof of Corollary 2.13. Here, the analysis is similar as in the case of continuous patterns, 
hence, we only emphasize the most important steps. 

Step 1. We fix a weak stationary solution (Uz, Vz, W%) with a null set J C [0, 1]. The 
Frechet derivative of the nonlinear mapping J 7 : Hz — > R-i defined by the mappings 
(fi, f2, i^) in (6.1) at the point (Uz, Vz, Wz) £ R-z has the form 



DT(Uz,Vz,W x ) V =Az(x 
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where Ax(x) = A(x) (see the matrix in (6.3)) if x G [0, 1] \ X and 

/ \ 
At(x) = if x £ I. 






This results immediately from the definition of the Frechet derivative. 
Step 2. Next, we study spectral properties of the linear operator 

ftp\ / \ 
C x V = + A 

\ v J v ° ° H 2 W 

in the Hilbert space U x (see (2.33)) with the domain D(C X ) = Lf (0, 1) x Lf (0, 1) x 
W 2,2 (0,1). First, we notice that the counterpart of Lemma 6.1 holds true and Ao (the 
positive eigenvalue of A\2 ) belongs to the continuous spectrum of (C x , D(Ci)). To prove 
this property of Ao, it suffices to follow the proof of Lemma 6.1. 

Step 3. A complex number A is an eigenvalue of the operator C x if and only if the 
following two conditions are satisfied: (i) the number A is not an eigenvalue of the matrix 
A\2', (ii) the boundary value problem 

-r J " + qx(;X)r 1 ^0, x G (0, 1) 
(6.22) 7 

7/(0) = 7/(1) = 

has a nontrivial solution. Here, we denote 

t xn f d X x) ~$ if xG [0,1]\X 

q x ( X , A) = < det(.4i2-A/) 

I — dg ~ ^ if x G X. 

This is the counterpart of Lemma 6.2 with an almost identical proof. 

Step 4- There exists a sequence {A n }„ G N of positive eigenvalues of the operator C x 
that satisfy A n — > Ao as n — > oo. This is the statement of Theorem 2.10 written for 
discontinuous patterns and its proof follows almost the same arguments. In particular, in 
Step 3 of the proof of Theorem 2.10, we should use the fact the the mapping A H- q x (-, A) 
is continuous in the L 2 (0, l)-norm (cf. Lemma 6.5). 

Step 5. Since, — C x is a sectorial operator, we complete the proof by [3, Thm. 5.1.3]. □ 



Appendix A. Kinetic system 

Here, we briefly review results on the large time dynamics of solutions to reaction- 
diffusion equations (2.1)-(2.3) supplemented with constant initial values, namely, we con- 
sider the following system of ordinary differential equations (so-called the kinetic system): 

. du ( av ,\ 

A.l _ 4 W 

at \u + v / 

, . n s dv _ _ 2 _ 

(A. 2) — = — dbV + u w — dv, 

i a on dw , _ _9 _ ,_ 

(A. 6) — — = —d g w — u w + dv + kq. 

(JjV 
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Steady states. It is easy to see that (A.1)-(A.3) has the trivial steady state (uo, Vo, Wq) = 
(0, 0, Ko/dg). Here, the right-hand side of equation (A.l) is satisfied in the limit sense, 
namely, when v \ and u \, 0. On the other hand, assuming that a > d c and 

d 2 c (d b + d) 



(A.4) k 2 Q > 6, where 6 = Ad g d b - 



[a - d c ) 



2 



system (A.1)-(A.3) has two positive equilibriums (u±, v±, w±), where 



_ a-d c _ _ d 2 c (d b + d) 1 _ k ± a/ k o ~ @ 

A.5 u± = — — v±, v± = j^- — , w± = — . 

d c (a — d c y w± 2d g 

Indeed, for m ^ 0, we obtain from equation (A.l) (with du/dt = 0) that u = v(a — d c )/d c . 
Substituting this expression to equation (A. 2) (with dv/dt = 0), we obtain 

d 2 c (d b + d) 1 



(A.6) 



{a — d r ) 2 w 



Finally, adding equations (A. 2) to equation (A. 3) (with dv/dt = dw/dt = 0) and using 
expression (A.6) we obtain the quadratic equation 

/ . _ % , 2 , d 2 (d b + d) 

A.7 d g w 2 - k w + d b ; v - = 0. 

(a - d c ) 2 

It is clear that equation (A.7) has two roots w± if Kq > B and one root w± = Ko/(2d g ) if 
k 2 = B, see (2.11)-(2.12). 

Boundedness of solutions. In the following, we denote by (u(t),v(t),w(t)) a solution 
of system (A.1)-(A.3) with a nonnegative initial datum. 

Proposition A.l. Every solution (u(t),v(t),w(t)) of system (A.1)-(A.3) corresponding 
to a nonnegative initial datum exists for all t > 0. Moreover, it is nonnegative, bounded 
for t > 0, and satisfies 

(A.8) limsupu(t) < — < liminf (v(t) + w(t)) < limsup (v(t)+w(t)) < — , 

i->oo [id c V t^oo n 

where v = max{d b , d g } > and \x = min{c4, d g } > 0. 

Proof. The global existence of nonnegative solutions is a consequence of Theorem 3.1. To 
show the relations in (A.8), it suffices to follow the proof of Theorem 2.1. □ 

Convergence to the trivial steady state. It is already proven (see Corollary 2.4 and 
Proposition 4.1) that the trivial steady state (0,0, K /d g ) is locally asymptotically stable 
as a solution to reaction diffusion system (2.1)-(2.3). In the case of the kinetic system, we 
can describe more precisely the convergence of solutions towards the trivial steady state. 

First, we consider the case a < d c , where, by Proposition 4.1, all solutions of the 
reaction-diffusion system converge to the trivial steady state. Now, we improve this 
result for system (A.1)-(A.3) 

Theorem A. 2. If a < d c , then every positive solution (u(t),v(t),w(t)) converges, as 
t — > oo, towards (0,0, Ko/d g ). Moreover, 

• if d c < d b + d, then v(t)/u(t) — >■ as t — > oo; 

• if d c > d b + d and a > d c — (d b + d) , then v(t)/u{t) — > d c (d b + d) / \d b + d + a — d c ) 
as t — > oo; 

• if d c > d b + d and a < d c — (d b + d), then v(t)/u(t) — > oo as t — > oo. 
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Next, we deal with the case a = d c , where the result is analogous. 

Theorem A. 3. If a = d c , then every positive solution (u(t),v(t),w(t)) converges, as 
t —> oo, towards (0, 0, n /d g ). Moreover, 

• if a < db + d, then v(t)/u(t) — > as t — > oo; 

• if a > db + d, then v(t)/u(t) — >• (a — (db + d)) / (db + d) as t — > oo. 

For a > d c , we see from Corollary 2.4 that the trivial steady state is locally asymp- 
totically stable. This convergence can be better described in the case of solutions of the 
kinetic system. Notice that if Kq < 0, there is no other nonnegative constant steady 
states of system (A.1)-(A.3). 

Theorem A. 4. Assume that a > d c and Kq < G. Then every positive solution 
(u(t),v(t),w(t)) converges, as t — >• oo, towards (0, 0, Ko/d g ). Moreover, 

• if d c < d b + d, then v(t)/u(t) — > as t — > oo; 

• if d c > d b + d, then v(t)/u(t) — > (d c — (d b + d)) / '(db + d + a — d c ) as t — )■ oo. 

If Kq > G, there are positive steady states defined in (A. 5). Letting an initial datum 
u(0) to be sufficiently small, we obtain that a solution converges to the trivial steady 
state, again. 

Theorem A. 5. Assume that a > d c and > G. Let (u(t),v(t),w(t)) be a solution to 
(A.1)-(A.3) corresponding to a positive initial condition (u(0),v(0),w(0)) . Ifu(0) satisfies 

_ a - d c K Q - a/kq - G 
M(0) < U+ " ~d c 2d b ' 

then this solution converges towards the trivial steady state. Here, u + is the constant 
defined in (A. 5). 

Since ideas of the proofs Theorems A.2-A.5 are more-or-less the same, we are going to 
sketch below the proof of Theorem A. 4, only. Introducing new variables X = v/u and 
Y = uw and writing, for simplicity, u instead of u, we obtain the following system 

... du ( aX , \ 

(A.10) ^ = -(d b + d + a-d c )X+-^— + Y, 

at 1 + X 

dY ( aX \ 

(A.ii) — = ( ~ rfc ) Y ~ dgY ~ U ^ Y + dXu2 + K ° u 

supplemented with positive initial conditions 

u(0) = u , X(0) = ^ = X , Y(0) = u(0)w(0) = Y . 

«(0) 

If u(t), v(t), and w(t) are positive of all t > 0, the functions X(t) and Y(t) are also 
positive. The following lemma plays an important role in the proof. 

Lemma A. 6. Let a > d c and Kq < G. If Xq < d c j (a — d c ) and Y < d c (db + d)/(a — d c ), 
then X(t) and Y(t) satisfy 

(A.12) 0<X(t)<— ^— and < Y(t) < ^ + , ^ for all < t < oo. 
v ' y ' a-d r w a - d r 



REACTION-DIFFUSION MODEL OF EARLY CARCINOGENESIS 



35 



Sketch of the proof of Lemma A. 6. Since this lemma can be shown by a method similar 
to that from the proof of Lemma 4.3, we just mention an important point, only. Here, we 
derive the following inequality 



(A.13) 

where 
(A. 14) 



Y(t) 



d c (d b + d) 



< -(Clg + Ul 



J t 



Y(t) 



d c (d b + d) 
a — d r 



R(u) 



R(u) 



d c d g (db + d) 
a — d c 



d c d b o , 

-u + k u. 



a 



If R(u) < 0, we immediately obtain Y(t) < d c (db + d)/(a — d c ) for all t > 0. However, it 
is easy to check that the assumption Kq < 9 implies R{u) < for all u > 0. □ 

Now, notice that, under the assumptions a > d c and Kq < O, the inequality X(t) < 
d c /(a — d c ) implies u t (t) < (see eq. (A. 9)). Consequently, if X < d c /(a — d c ) and 
Y < d c (db + d)/(a — d c ), Lemma A. 6 provides an estimate of a solution (u(t), v(t), w(t)) 
for t > 0, which allows us to show its convergence towards the trivial steady state (see 
the proof of Theorem 2.2). On the other hand, if either X Q > d c /(a — d c ) or Y Q > 
d c {db + d)/(a — d c ), one can show, by simple geometric arguments involving equations 
(A.9)-(A.10) and the boundedness of solutions to system (A.1)-(A.3), that there exists 
T > such that X{T) < d c /(a — d c ) and Y(T) < d c (db + d)/(a — d c ). It is easy to see 
that every solution of (A.9)-(A.ll) stays in the rectangle (A. 12). 

In order to show the convergence rate of u and v which are stated in Theorem A. 4, we 
consider the steady states of system (A.9)-(A.ll). It is clear that (u,X, Y) = (0,0,0) is 
an equilibrium. Moreover, noting that Y(t) — > as t — > oo when u(t) — > as t — > oo, we 
have the following nonnegative steady state of system (A.9)-(A.ll), 



(A.15) 



(u,X,Y)= 



d c — (db + d) 



db + d 



d r 







if d c > db + d. 



In the following, we denote by E\ the right-hand side of (A.15). We study stability of E\ 
by analyzing eigenvalues of the corresponding Jacobian matrix at the equilibrium. 

Lemma A. 7. If d c < db + d, then the equilibrium (0,0,0) is asymptotically stable, while 
it is unstable if d c > d b + d. For d c > db + d, the stationary solution E\ is asymptotically 
stable. 

Proof. The Jacobian matrix J of the nonlinear mapping defined by the right-hand side of 
system (A.9)-(A.ll) is of the form 



J 



V 



aX 
1+X 



— d r 





-2uY + 2dXu + k 



-(d b + d + a- d c ) + (1+ a x)2 



aY 



+ du 1 



aX 
1+X 




1 

d c -d g - u 2 



Hence, at the stationary solution (0,0,0), it becomes 



J, 



(0,0,0) 






-(db + d + a — d c ) 
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and it is clear that all eigenvalues of J(o,o,o) have negative real parts. Similarly, letting 
Je 1 to be the corresponding Jacobian matrix at the steady states E±, we obtain 



J 










(d b +d+a— d c ){d b +d— d c ) 

-(d b + d + d g ) 



a 





It follows from inequality d c > db + d that all eigenvalues of Je 1 are negative. 
Now, the assertions in Theorem A. 4 follow from Lemmas A. 6 and A. 7. 



□ 



Remark A. 8. For > O, we have R(u) < (cf. (A. 14)) provided u < u + . Therefore, 
under the condition u(0) < u + , we obtain the same assertion as in Lemma A. 6, now 
however, under the assumption a > d c and Kq > 0. Here, it suffices to use the fact 
that du/dt < whenever X(t) < d c /(a — d c ). Repeating the reasoning from the proof of 
Theorem 2.2, we can show Theorem A. 5. 

Stability of positive steady states. Let (u±,v±,w±) be the positive steady states of 
system (A.1)-(A.3) given by formulas (A. 5). To study their stability, we consider again 
system (A.9)-(A.ll), which for a > d c and kq > 0, has two positive steady states 

d c d c (d b + d) 



(u,X,Y) ± = (u ± , 



a — d c a — d c 



Theorem A. 9. The vector (u,X,Y)_ is an asymptotically stable stationary solution of 
system (A.9)-(A.11) ; and (u,X,Y) + is unstable. 

Proof. First we consider Jacobian matrix J + of the right-hand side of (A.9)-(A.ll) at 
(u,X,Y) + : 



( 



J- 








(q-rfc) 2 - 



d b + d + 



d c (a—d c )(d b +a 



U + 
dc(a—d c ) 



+ du\ 





1 



-d n - 2 



-9 K ) 

Hence, every eigenvalue A of J + is determined by its characteristic equation F(\) = 0, 
where 



F(\) 



d b + d+ d q + 



dJa - d r 



A 



(A.16) 



(d g + u + ) I d b + d + 



d c (a - d c ) 



d c (a 



d r 



(d b + d) + du 



-2 



llfi 



0- 



d r 



Since F(0) > and F(+oo) = — oo, there exists at least one positive real eigenvalue. This 
implies that (u, X, Y) + is unstable. 

Next, we investigate Jacobian matrix J_ at (u, X, F)_, 



/ 








(a-d c ) 2 - 



d b + d + 



d c (a—d c ) 



^0 



d c (a-d c )(d b +d) 
a 



+ du 



v, 2 





1 



-d - 2 
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Here, calculations are more involved to show stability of (u,X,Y)_. We consider again 
the characteristic equation G(X) = 0, where 



G(A) 



-A a 



a 2 



d b + d + d„ + 



dJa - d r 



(A.17) 



(dg 



dh + d 



d c (a - d c ) 



A 2 

d c (a — d c 
a 



— \ K 



(a-d c f 



(7 



Now, for simplicity, we set G(X) 



-A 3 - AA 2 - BX - ^4 



■(d b + d) + dut 



9 (a - 4) V 



Since 



A 2 — 3B > 0, the equation G"(A) = has two negative roots 



Ai 



A + vM 2 - W 



A, 



A - vM 2 - 3B 



-3 ' -3 

This implies that G(X) = has a local minimum and a local maximum in the half 
plane {A : A < 0}. Moreover, we obtain that G(Ai) < by virtue of the inequality 
A 2 -3B > 0. Therefore, noting that G(Q) < 0, G(-oo) = +oo and G(+oo) = -oo, we 
have two possibilities: 

either (i) G(A 2 ) > or (zi) G(A 2 ) < 0. 

Here, G(A 2 ) is given by 



G(A 2 



A 2 



A + 2VA 2 - 3B 



n 2 -Q 



'a — d r ) 2 



If (i) occurs, then the equation G(X) = has three negative real roots, which implies 
immediately that (u,X, y)_ is asymptotically stable. Here, we have to notice that the 
case (i) occurs if is sufficiently close to 6. 

Next, suppose G(X) satisfies (ii), namely, G(A 2 ) < 0. Then, the equation G(A) = 
has three roots: one is real and negative, we denote it by — Ao, and others are complex 
numbers fi and ft. Therefore, using the decomposition G(A) = — (A + A )(A — /S)(X — ft), 
we obtain from (A.17) that 



A = A - 2Re(/i), B 
Now, we consider the following quantity 

AB - 



|/4 1 — A 2Re (n). 



(A.18) 



G{-A) 



e (a - 4) V 



Differentiating the right-hand side of (A.18) with respect to k,q, we obtain 
d 



(A.19) 



dno 



AB 



9 (a - 4) V 



> for all nl > 0. 



Since G{— A) > at = 6, inequality (A.19) implies that G(— A) > for all > 
©. Under the assumption (ii), we notice that the inequality G(X) > is satisfied only 
when A < — Ao- Therefore, we obtain the inequality —A < — Ao, which is equivalent to 
A — 2Re (/i) > A . This implies Re (fi) < 0. Consequently, if the case (ii) occurs, then real 
parts of all eigenvalues of J_ are negative, which implies that (u,X, F)__ is asymptotically 
stable. □ 
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Since X = v/u and Y = uw, the steady states (u,X,Y)_ and (u,X,Y) + of system 
(A.9)-(A.ll) correspond to the steady states (u_,v^, wJ) and (u + ,v + ,w + ) of (A.l)- 
(A.3), respectively. Here, one should recall the expressions from (A. 5) to see that u± = 
d c (db + d)/((a — d c )w±), hence, 

a — d c k + \/ kI — ® a — d c k — \/ — 

u - = ~i h and u + = ~a h • 

d c 2d b d c 2db 

Let us also notice that U- > u + . 

Thus, we have proved that stationary solutions of the original kinetic system (A.l)- 
(A.3) have the following stability properties. 

Corollary A. 10. The vector is an asymptotically stable stationary solution 

of system (A.1)-(A.3), while (u + ,v + ,w + ) is unstable. 
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